Add DAEProblem integration

[ChrisRackauckas]

using OrdinaryDiffEq
using DiffEqGPU

function f1(out,du,u,p,t)
    out[1] = p[1]*(u[2]-u[1]) - du[1]
    out[2] = u[1]*(p[2]-u[3]) - u[2] - du[2]
    out[3] = u[1]*u[2] - p[3]*u[3] - du[3]
end

u₀= Float32[1.0;0.0;0.0]
du₀ = Float32[0.5;0.5;0.5]
tspan = (0.0f0,100.0f0)
p = [10.0f0,28.0f0,8/3f0]
differential_vars = ones(3)

prob = DAEProblem(f1,du₀,u₀,tspan,p,differential_vars = differential_vars)
prob_func = (prob,i,repeat) -> remake(prob,p=rand(Float32,3).*p)
monteprob = EnsembleProblem(prob, prob_func = prob_func, safetycopy=false)

@time sol = solve(monteprob,DFBDF(),EnsembleGPUArray(),trajectories=1_000,saveat=1.0f0)