Debugging differential equation solver issues almost always boils down to doing the same thing, so this is a summary to help you out. For more information on specific issues, check out the FAQ section of the DifferentialEquations.jl documentation which highlights common issues and questions:

https://diffeq.sciml.ai/dev/basics/faq/

## How do I debug why the differential equation solver is diverging?

dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.

NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.

Instability detected. Aborting

There, now this post will show up in a search. What do you do when you see these kinds of problems?

If you find yourself in differential equation solver trouble, the first thing to do is lower the tolerance. In many cases, decreasing the tolerance improves the stability, and so trying something like `abstol=1e-10,reltol=1e-10`

can tell you if itâ€™s a tolerance issue. If thatâ€™s not working, try switching to a more stable solver, something for stiff equations like `TRBDF2()`

, `KenCarp4()`

, or `QNDF()`

.

If none of those things are working, the problem is probably your model. If you want to blame Julia solvers, then thereâ€™s a good way to prove it: try a non-Julia solver. Itâ€™s as simple as changing `solve(prob,Tsit5())`

to `solve(prob,CVODE_BDF())`

and now youâ€™re using the classic Sundials C/C++ library. Thatâ€™s not all, thereâ€™s a whole list:

- Sundials.jl, a wrapper for the C/C++ SUNDIALS library though
`CVODE_Adams`

,`CVODE_BDF`

,`IDA`

, and`ARKODE`

. - ODEInterfaceDiffEq.jl, a wrapper for the classic Hairer Fortran codes like
`dorpi5`

,`dop853`

,`radau`

,`rodas`

, etc. - LSODA.jl, a wrapper for the classic
`lsoda`

algorithm. - MATLABDiffEq.jl, a wrapper for the MATLAB ODE solvers
`ode45`

,`ode15s`

, etc. - SciPyDiffEq.jl, a wrapper for SciPyâ€™s
`odeint`

(LSODA) and other methods (LSODE, etc.). - deSolveDiffEq.jl, a wrapper for the commonly used R library.

*Note: these solver packages are not included by default due to their extra dependencies, so youâ€™lld need to do things like ]add Sundials; using Sundials before using the Sundials solvers*

If you do this at lower tolerances and every solver failsâ€¦ Iâ€™ll just put it in meme form.

So yes, if youâ€™re model is failing with every major solver, including all of the major ones called into unmodified from C/C++ and Fortran, then the problem isnâ€™t the solver itâ€™s your model. I know you checked it, I know you think you debugged it, but the likelihood that every single solver ever created in all languages is incorrect instead of the code you wrote a few hours ago is very slim.

So what can it be? Hereâ€™s a list of common issues:

- Double check for sign issues. And then check again. If you have some solution piece thatâ€™s blowing up to infinity, why? Which term in its derivative gets really large? Is it supposed to be large? Why did it get large? Almost all issues boil down to just checking term by term in this fashion. Double check your model.
- Double check your assumptions on the model. If the derivative goes to zero as
`u`

goes to zero, that means that`u`

can never go negative, right? Wrong!`u' = -sqrt(u)`

hits zero in finite time. Just because you believe that the system youâ€™re modeling has a property (like being positive) doesnâ€™t mean that the true solution to your model actually has this property. Try to prove the properties or, try and look at the solution as it begins to violate the property, look at the values in the derivative, and understand which terms are causing it. - Double check if youâ€™re doing anything that violates assumptions of being an ODE. As an ODE, the right-hand side
`f`

function should always give you the same result:`u' = f(u,p,t)`

needs to be uniquely defined. These issues are just fundamental to the mathematics: if you do these things in`f`

, then`f`

no longer defines an ODE so of course it cannot be solved!- If you put randomness into
`f`

, then the adaptivity will think your ODE is being solved at a high error because the derivative keeps changing, and therefore it will fail and hit`dtmin`

trying to reduce the randomness to zero (if you do need randomness, use an SDE or RODE solver). - If your
`f`

function is modifying`u`

, then calling`f`

with different stepsizes is not deterministic and the solver is likely to fail. If you need to do this, you should be using a callback. - If your
`f`

function is caching values from a previous step, just think about what this means. If you change`dt`

, then youâ€™re changing`f`

. In that sense,`u'`

is no longer defined since itâ€™s now dependent on how itâ€™s being solved! Even worse, adaptive ODE solvers do not always move forwards in time: sometimes they try`t + dt1`

before trying`t + dt2`

and choosing what to do. So if youâ€™re assuming that the cached values are from the last step, thatâ€™s not actually the case: those values can be coming from a fake (too incorrect according to error estimates) future!

- If you put randomness into

If you â€śknowâ€ť your model canâ€™t have an issue because itâ€™s translated from MATLAB/Python/R, then try using MATLABDiffEq.jl/SciPyDiffEq.jl/deSolveDiffEq.jlâ€™s wrapper solvers, which call the solvers of those respective languages/packages. If your differential equation is the same, then it should give exactly the same results, down to the step. Almost certainly if the Julia solvers are diverging these will diverge as well, and you can use this to help narrow down translation problems.

## What about performance issues?

For understanding how to solve differential equations with the top-notch performance, follow the tutorials which mention how to get high performance:

https://diffeq.sciml.ai/dev/tutorials/advanced_ode_example/

https://tutorials.sciml.ai/html/introduction/03-optimizing_diffeq_code.html

https://mitmath.github.io/18337/lecture2/optimizing

If youâ€™ve gone through these tutorials but still have performance issues, or just have clarifying questions to ask, please feel free to ask! But do check out these tutorials since they cover most of what people need!

## Benevolent onlookers: how to help

Thereâ€™s a large volume of differential equation questions pouring in because thereâ€™s a lot of adoption in undergraduate and graduate courses. This is a good thing! But itâ€™s also time consuming and >99% of questions come down to the issues listed above and have a formulaic solution. If youâ€™ve read this post, you now have enough information to lead someone through debugging model issues so please feel free to share this and walk someone through the step by step debugging process. If it turns out to be one of the few cases that need extra attention, then point people towards opening an issue.