Solving same DiffEq repeatedly but with different initial conditions

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]))

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)]

# How do you get rid of this whole loop?
tspan = (0.0,abs(minimum(
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
# do some more processing on the saved result for all x

You can use the integrator interface with reinit.

Or you can use the ensemble interface to auto-parallelize it.

This includes automated GPU parallelism with DiffEqGPU.jl


Thanks Chris, that resolves it for my use case.

That’s wonderful, I’ll definitely move to this at some point if my system gets large enough.