What is the best way to solve so many decoupled ODEs? Possibly in parallel.

Should I use Modeling Toolkit or DifferentialEquations?

The equations are of the type: u’ = -a^2*k^2*u

where a is the same for all, then I have a vector with the values of k and a vector with the initial conditions u0.

Hi, welcome to Julia community. Although it is, of course, completely up to you how you pose your question/request, I dare to suggest that you first implement whatever solution and share the code here (this is also the recommended procedure here). My experience is that suggestions and advice come better focused then. You (and everybody else) can also immediately check if the advice leads to an improvement with respect to the initial “non-expert” solution.

The differential equation looks linear and thus has a closed-form solution. You could simulate all equations at once with some vectorized calls to the exponential function

```
a = 1.0
n = 5
K = 1:n
u0 = rand(n)
t = range(0.0,1.0,100)
u = [exp.(-a^2*k^2*t) for k in K].*u0
using Plots
plot(t,u)
```

These are my attempts:

```
function solution1(a, ys, ks, t_start, t_max)
f(u, p, t) = @. -a^2*p^2*u
prob = ODEProblem(f, ys,(t_start, t_max), ks)
solve(prob, Tsit5())[end]
end
function solution2(a, ys, ks, t_start, t_max)
sol = Array{typeof(ys[1])}(undef, length(ys))
for i in eachindex(sol)
@parameters t α κ
@variables u(t)
∂t = Differential(t)
eq = [∂t(u) ~ -α^2*κ^2*u]
@named sys = ODESystem(eq)
u0 = [u => ys[i]]
p = [α => a, κ => ks[i]]
prob = ODEProblem(sys, u0, (t_start, t_max), p)
sol[i] = solve(prob, Tsit5())[end][1]
end
sol
end
```

solution2 is much faster than solution1. For example if ks and ys are Complex64 vectors (where length(ks) = length(ys) = 700 ) i get:

30.545022 seconds (3.19 M allocations: 33.188 GiB, 56.47% gc time) for solution1

10.310859 seconds (88.86 M allocations: 7.452 GiB, 20.72% gc time) for solution2

If it’s just a linear ODE, follow the advice and just use `exp`

.