I have a set of coupled ODEs that I’d like to solve, which is simple enough using DifferentialEquations.jl. My problem is that I need to solve the same ODE repeatedly for many different initial conditions parametrized by an array x.
Currently I just do this in a loop. I’m sure there’s a smarter way but I haven’t quite found it in the documentation yet. This is very inefficient of course since I need to keep setting up the problem again every time.
Could someone give me some pointers?
Here is a minimal non-working example of what I’m doing:
using Elliptic.Jacobi
function dn_ab!(du,u,p,t)
du[1] .= +im*p[1].*u[1] .+ im.*u[2].*dn.(t,real(p[2]))
du[2] .= -im*p[1].*u[2] .+ im.*u[1].*dn.(t,real(p[2]))
end
function ab_dn_t0(x)
A = +χ .+ 0.5.*Ω*λ.*(x - xs) .- π/4
B = -χ .+ 0.5.*Ω*λ.*(x - xs).- π/4
return [2*im*sin(A); 2*cos(B)]
end
# How do you get rid of this whole loop?
tspan = (0.0,abs(minimum(c.box.t)))
for i in eachindex(x)
u0 = ab_dn_t0(x[i])
prob = ODEProblem(dn_ab!, u0, tspan, [λ, m])
sol = solve(prob, Tsit5(), saveat=t) # t is an array
# do some processing on sol and save the result
end
# do some more processing on the saved result for all x