# Solve many first-order ODEs

**URL:** <https://discourse.julialang.org/t/solve-many-first-order-odes/85060>\
**Category:** Modelling & Simulations\
**Created:** [July 31, 2022, 1:48pm UTC](https://discourse.julialang.org/t/solve-many-first-order-odes/85060 "2022-07-31T13:48:20Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![qwerty](https://avatars.discourse-cdn.com/v4/letter/q/4491bb/32.png) [@qwerty](https://discourse.julialang.org/u/qwerty)\
**Post date:** [July 31, 2022, 1:48pm UTC](https://discourse.julialang.org/t/solve-many-first-order-odes/85060/1 "2022-07-31T13:48:20Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [July 31, 2022, 2:02pm UTC](https://discourse.julialang.org/t/solve-many-first-order-odes/85060/2 "2022-07-31T14:02:47Z")

</div>

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](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757) 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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 31, 2022, 2:11pm UTC](https://discourse.julialang.org/t/solve-many-first-order-odes/85060/3 "2022-07-31T14:11:36Z")

</div>

[https://diffeq.sciml.ai/stable/features/ensemble/](https://diffeq.sciml.ai/stable/features/ensemble/)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 31, 2022, 2:19pm UTC](https://discourse.julialang.org/t/solve-many-first-order-odes/85060/4 "2022-07-31T14:19:04Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [July 31, 2022, 3:57pm UTC](https://discourse.julialang.org/t/solve-many-first-order-odes/85060/5 "2022-07-31T15:57:10Z")

</div>

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![qwerty](https://avatars.discourse-cdn.com/v4/letter/q/4491bb/32.png) [@qwerty](https://discourse.julialang.org/u/qwerty)\
**Post date:** [July 31, 2022, 4:37pm UTC](https://discourse.julialang.org/t/solve-many-first-order-odes/85060/6 "2022-07-31T16:37:10Z")

</div>

These are my attempts:

```julia
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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 31, 2022, 6:25pm UTC](https://discourse.julialang.org/t/solve-many-first-order-odes/85060/7 "2022-07-31T18:25:21Z")

</div>

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