# Choosing suitable options for solver ODE

**URL:** <https://discourse.julialang.org/t/choosing-suitable-options-for-solver-ode/77843>\
**Category:** Numerics\
**Tags:** question, diffeq, ode\
**Created:** [March 14, 2022, 3:07am UTC](https://discourse.julialang.org/t/choosing-suitable-options-for-solver-ode/77843 "2022-03-14T03:07:08Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Camilo\_Chaparro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/camilo_chaparro/32/33597_2.png) [@Camilo\_Chaparro](https://discourse.julialang.org/u/Camilo_Chaparro)\
**Post date:** [March 14, 2022, 3:07am UTC](https://discourse.julialang.org/t/choosing-suitable-options-for-solver-ode/77843/1 "2022-03-14T03:07:08Z")

</div>

Hello everyone, I want about the most suitable ODE solver for the following problem:

```julia

using DifferentialEquations
using Plots; plotly() 

#Parameters

omega_m = 1.6e6         
N = 1
Delta_mw = -0.7           
g0 = 3.6e3#/omega_m 
g_eff = (g0/sqrt(N))
Gamma = 60000#/omega_m
Gamma_tilde = 0.0       
Omega_mw = 0.7           
Q = 2*(10^5)
gamma_m = (omega_m/Q) #/ omega_m
N_th = 6.05
Angle_th = 0;
tau = (Gamma/(g_eff^2));
#omega_m = 1 

function mean_field(du,u,p,t)
    
    
    #p[1]=Delta_mw
    #p[2]=Omega_mw
      
    Sx =u[1]
    Sy =u[2]
    Sz =u[3]
    alpha_x=u[4]
    alpha_y=u[5]
        
    du[1] = (Delta_mw-4*g_eff*u[4])*u[2] - (Gamma+Gamma_tilde)*u[1]
    du[2] = -(Delta_mw-4*g_eff*u[4])*u[1] -Omega_mw*u[3] - (Gamma+Gamma_tilde)*u[2]
    du[3] = Omega_mw*u[2] - 2*Gamma*(u[3]+N)
    du[4] = omega_m*u[5] - (gamma_m /2)*u[4]
    du[5] = -omega_m*u[4] -g_eff*u[3] - (gamma_m /2)*u[5]

    
    
    
    #du[1] = (p[1]-4*g_eff*alpha_x)*Sy - (Gamma+Gamma_tilde)*Sx
    #du[2] = -(p[1]-4*g_eff*alpha_x)*Sx -p[2]*Sz - (Gamma+Gamma_tilde)*Sy
    #du[3] = p[2]*Sy - 2*Gamma*(Sz+N)
    #du[4] = omega_m*alpha_y - (gamma_m /2)*alpha_x
    #du[5] = -omega_m*alpha_x - g_eff*Sz - (gamma_m /2)*alpha_y
    
    end  

function Result_lim(Array_sol)
    
    print((Array_sol[2][4])^2+(Array_sol[2][5])^2," ")
end

#Initial state

u0 = [0;0;-N;sqrt(N_th);0];
time= 100*tau

function Quanta_fin(Delta_mw,Omega_mw)
    
    parameters=[Delta_mw ;Omega_mw]
    tspan = (0.0,time)
    prob = ODEProblem(mean_field,u0,tspan,parameters) #,parameters)
    
    Val1=1e-6
    Val2=1e-6
    Iter=1e14
    dt_size=1e-7
 
    sol1 = solve(prob,RK4(),dt=dt_size,save_everystep=false,reltol=Val1,abstol=Val2,MaxIters=Iter)
    sol2 = solve(prob,alg_hints=[:stiff],dt=dt_size,reltol=Val1,abstol=Val2,save_everystep=false,MaxIters=Iter)
    sol3 = solve(prob,Rodas4() ,dt=dt_size,reltol=Val1,abstol=Val2,save_everystep=false,MaxIters=Iter)
    sol4 = solve(prob,EulerHeun() ,dt=dt_size,reltol=Val1,abstol=Val2,save_everystep=false,MaxIters=Iter)
    sol5 = solve(prob,OwrenZen5() ,dt=dt_size,reltol=Val1,abstol=Val2,save_everystep=false,MaxIters=Iter)
    sol6 = solve(prob,TanYam7() ,dt=dt_size,reltol=Val1,abstol=Val2,save_everystep=false,MaxIters=Iter)
    sol7 = solve(prob,alg_hints=[:nonstiff] ,dt=dt_size,reltol=Val1,abstol=Val2,save_everystep=false,MaxIters=Iter)

    
    Array_sol1=sol1.u
    Array_sol2=sol2.u
    Array_sol3=sol3.u
    Array_sol4=sol4.u
    Array_sol5=sol5.u
    Array_sol6=sol6.u
    Array_sol7=sol7.u

    Result_lim(Array_sol1)
    Result_lim(Array_sol2)
    Result_lim(Array_sol3)
    Result_lim(Array_sol4)
    Result_lim(Array_sol5)
    Result_lim(Array_sol6)
    Result_lim(Array_sol7)

    
    end

Quanta_fin(-Delta_mw,Omega_mw)

```

I tried to use multiple solvers and even I decreased the step size and increased the number of iterations and as you can see the results are different (Rodas4() and alg\_hints=[:stiff] have similar values), the main problem relies long time needed to achieve the steady state dynamics. That is why I want to ask how can I change the options of the solver to obtain a reliable solution and predict properly the steady state of the system. Thanks.

---

<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:** [March 14, 2022, 1:19pm UTC](https://discourse.julialang.org/t/choosing-suitable-options-for-solver-ode/77843/2 "2022-03-14T13:19:24Z")

</div>

If you have a lower error tolerance for the solution, then decrease the error tolerance. And this would run hundreds of times faster if you didn’t use non-constant globals.

---

<div class="post-metadata">

**Author:** ![Brad\_Carman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brad_carman/32/17631_2.png) [@Brad\_Carman](https://discourse.julialang.org/u/Brad_Carman)\
**Post date:** [August 6, 2025, 1:28pm UTC](https://discourse.julialang.org/t/choosing-suitable-options-for-solver-ode/77843/3 "2025-08-06T13:28:41Z")

</div>

I added this example to my new analysis package: [bradcarman/ModelingToolkitTolerances.jl: A package to help determine `abstol` and `reltol` when using ModelingToolkit.jl and DifferentialEquations.jl](https://github.com/bradcarman/ModelingToolkitTolerances.jl)

See [ModelingToolkitTolerances.jl/examples/777843.jl at main · bradcarman/ModelingToolkitTolerances.jl](https://github.com/bradcarman/ModelingToolkitTolerances.jl/blob/main/examples/777843.jl) for a detailed analysis of a few solvers. Finding the right solver still requires a lot of knowledge, but with this package it helps to quickly explore solver choices along with a variation in tolerance settings to more quickly locate the best option. Here I find that `Rodas4` does the best, with a max residual of 0.176 at abstol=1e-9 and reltol=1e-12.

Here is what the Rodas4 analysis looks like…

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/b/1b34ab15396154fe904ae10315209267222488a7.png)

```julia-auto
┌────────┬────────────────┬─────────────────┬─────────────────┬──────────────────┐
│ abstol │ reltol = 0.001 │ reltol = 1.0e-6 │ reltol = 1.0e-9 │ reltol = 1.0e-12 │
├────────┼────────────────┼─────────────────┼─────────────────┼──────────────────┤
│ 0.001 │ 16310.4 │ 3496.79 │ 1827.51 │ 1872.54 │
│ 1.0e-6 │ 12888.2 │ 65.512 │ 28.823 │ 26.6936 │
│ 1.0e-9 │ 16463.4 │ 82.4328 │ 0.440874 │ 0.175769 │
└────────┴────────────────┴─────────────────┴─────────────────┴──────────────────┘

```

Here is what the result looks like with the optimal settings…

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/9/595af5aee26b83d8ef6dfaa8e455276498905cec.png)