BatchIntegralFunction not appearing to work as intended

I am trying to implement a batch integral function using the new BatchIntegralFunction rather than then deprecated nout/batch keywords. However, when I do, the in-place vector to fill place the answer in is the incorrect size therefore it returns a broadcasting error. Please can you take a look at the MWE below and see where I went wrong. Thank you very much

function f(y, u, p)
    Threads.@threads for i in 1:length(u)
        y[1, i] = sin(u[i])*p[1]
        y[2, i] = cos(u[i])*p[2]
    end
end
p=[1,2]
int=BatchIntegralFunction(f,zeros(2),max_batch=2)
prob = IntegralProblem(int, 1, 3,p)#; nout = 2,nbatch=2
sol = solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)
sol.u

Open an issue.

Hi there, the prototype you passed to the BatchIntegralFunction is missing a batching dimension. Specifically, the prototype should have the same type as y, and the same shape along all dimensions, although the batching dimension is unspecified and can have length zero. The following works for me:

function f(y, u, p)
    Threads.@threads for i in 1:length(u)
        y[1, i] = sin(u[i])*p[1]
        y[2, i] = cos(u[i])*p[2]
    end
end
p=[1,2]
int=BatchIntegralFunction(f,zeros(2, 0),max_batch=2)
prob = IntegralProblem(int, 1, 3,p)#; nout = 2,nbatch=2
sol = solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)
sol.u

This is a point which is lacking in the Integrals.jl documentation so I am opening a pr.