OrdinaryDiffEq.jl: the default ODE solver of the SciML stack

A component package that most people meet as the whole ecosystem

OrdinaryDiffEq.jl is a component package in the DifferentialEquations ecosystem, and it holds the ordinary differential equation solvers and utilities. The readme is careful about what that means for a newcomer: the package is completely independent and usable on its own, but users interested in this functionality should check out DifferentialEquations.jl.

That distinction is the first thing to sort out before writing any code. DifferentialEquations.jl is the convenience layer that pulls in a large collection of solvers and helper packages. OrdinaryDiffEq.jl is the part that actually solves ODEs, and it also covers differential-algebraic equations, neural ordinary differential equations and the wider scientific machine learning area, which is the subject of the package description.

Installation is the standard Julia package manager call:

import Pkg;
Pkg.add("OrdinaryDiffEq");

The readme also notes the requirement that the user passes an OrdinaryDiffEq.jl algorithm to solve, which is a small but important detail. There is no default solver silently chosen for you, so the algorithm is part of your code and part of your decision.

The project sits in a mature organisational structure. There is a SciML common interface that the package implements, a ColPrac contributor guide badge, a SciML code style guide, CI through GitHub Actions plus a Buildkite pipeline, codecov coverage, a Zulip chat channel for the SciML bridged community, a CITATION.bib, a NEWS.md, a benchmark directory and a test directory. The documentation is not in the repository readme, it lives at docs.sciml.ai under DiffEqDocs.

Out-of-place for one equation, in-place for a system

The API example in the readme is worth reading closely, because the choice of function form is the first performance decision in Julia ODE work.

The out-of-place form returns a new value:

using OrdinaryDiffEq
f(u, p, t) = 1.01 * u
u0 = 1 / 2
tspan = (0.0, 1.0)
prob = ODEProblem(f, u0, tspan)
sol = solve(prob, Tsit5(), reltol = 1.0e-8, abstol = 1.0e-8)

The in-place form, described in the readme as more efficient for systems of equations, writes into a preallocated buffer:

using OrdinaryDiffEq
function lorenz!(du, u, p, t)
    du[1] = 10.0 * (u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
    return
end
u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz!, u0, tspan)
sol = solve(prob, Tsit5())

The pattern is consistent across the package. You build an ODEProblem from a function, an initial condition and a time span, then pass it to solve along with the algorithm and whatever tolerances you want. The adaptive method in these examples is Tsit5, and the tolerance arguments reltol and abstol are set explicitly in the first example.

For a model of a fixed size, the readme points at a third form that is described as very fast static array versions specifically compiled to the size of your model, using StaticArrays:

using OrdinaryDiffEq
using StaticArrays: SA
function lorenz(u, p, t)
    return SA[10.0 * (u[2] - u[1]), u[1] * (28.0 - u[3]) - u[2], u[1] * u[2] - (8 / 3) * u[3]]
end
u0 = SA[1.0; 0.0; 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz, u0, tspan)
sol = solve(prob, Tsit5())

The three forms are the same mathematics with different allocation profiles, and choosing between them is a per-model decision rather than a style preference.

Refined problems, symplectic methods and the SecondOrderODEProblem

Beyond plain first-order ODEs, the package handles refined ODE types, and the readme points at the DiffEqDocs pages on ODE types for the full list.

The example given is the harmonic oscillator, written as a second-order problem with the acceleration function taking the derivative buffer, the velocity, the position, the parameters and the time:

using OrdinaryDiffEq
function harmonic_oscillator!(dv, v, u, p, t)
    ω = p[1]
    dv[1] = -ω^2 * u[1]
    return
end
ω = 2.0  # angular frequency
initial_position = [1.0]
initial_velocity = [0.0]
tspan = (0.0, 10.0)
prob = SecondOrderODEProblem(harmonic_oscillator!, initial_velocity, initial_position, tspan, [ω])
sol = solve(prob, VelocityVerlet(), dt = 1 / 100)

The mechanism is a second-order reformulation. The equation is written in the readme as a second-order form and then shown to be equivalent to a first-order system, and SecondOrderODEProblem takes the initial velocity and the initial position in that order, with the parameters last.

The solver choice is the point of this example. VelocityVerlet is a symplectic method, and the readme explains why that matters: for complex dynamical systems such as the Henon-Heiles potential, symplectic integrators preserve the structure of Hamiltonian dynamics. The corresponding example uses the KahanLi8 method:

function HH_acceleration!(dv, v, u, p, t)
    x, y = u
    dv[1] = dx = -x - 2 * x * y
    dv[2] = dy = y^2 - y - x^2
    return
end
initial_positions = [0.0, 0.1]
initial_velocities = [0.5, 0.0]
prob = SecondOrderODEProblem(HH_acceleration!, initial_velocities, initial_positions, tspan)
sol2 = solve(prob, KahanLi8(), dt = 1 / 10);

Note that these symplectic solves pass a fixed dt rather than relying on adaptive stepping, which is inherent to the method class rather than an oversight in the example.

The readme also names IMEX and semi-linear ODEs, which is where exponential integrators come in, for problems where a stiff component and a non-stiff component are separated explicitly.

The full solver list is not in the readme. It is split across three documentation pages, one for ODE solvers, one for dynamical ODE solvers and one for split ODE solvers, which is a fair hint at how many algorithms are available.

The v7 upgrade touches everything underneath

The one piece of news in the readme that a user has to act on is the v7 breaking change.

OrdinaryDiffEq v7 bumps to SciMLBase v3 and RecursiveArrayTools v4, with breaking changes across all sublibraries. The readme does not attempt to summarise them, it points at NEWS.md in the repository as the full migration guide.

The phrase across all sublibraries is the important part. This is not a single package with a renamed function, it is a coordinated major version bump across the packages the solver is built on, which means the blast radius of the upgrade is wider than the diff in OrdinaryDiffEq itself. If your project pins versions across the SciML stack, this is the kind of upgrade that wants a branch and a test run rather than a manifest edit.

The dependency versions of the underlying stack matter for another reason. The choice of SciMLBase determines the problem and solution types, and the choice of RecursiveArrayTools determines how the solution array is stored and indexed, so both are part of your public surface whether you touch them or not.

The repository layout tells you where things live. There is a src directory for the main package, a lib directory that holds the sublibraries, a docs directory, a test directory, a benchmark directory for performance tracking, and a Project.toml at the root describing the package. The .typos.toml file is a small detail that indicates a check for spelling in the source and docs, which matters for a package whose identifiers are often mathematical notation.

The licence is worth flagging, and it is not the permissive one most people assume. The repository carries a LICENSE.md, and the licence is recorded as a custom one that GitHub cannot classify, so you need to read that file yourself before shipping it in a commercial or closed source product.

What this means if you are coming from another numerical stack

The reasonable alternative for a Julia user is to stay with DifferentialEquations.jl, which bundles the solvers together with the surrounding tooling rather than importing one component package. If you do not have a specific reason to control the dependency set, the readme itself points that way, and the practical difference is one import line versus a slightly smaller dependency tree.

The alternative for a user in another language is a solver from that language's ecosystem, and the trade is not only the solver. The advantage of OrdinaryDiffEq is not a single algorithm but the collection: adaptive methods, symplectic integrators for Hamiltonian systems, exponential integrators for stiff problems, split methods, and event handling for problems with discontinuities, all behind one common interface. Moving to another ecosystem usually means moving your model, your solver and your analysis code.

The other reason to be in Julia specifically is the neural ODE and scientific machine learning angle, which the package description names explicitly. If your project is differentiable-programming work where the ODE solver sits inside a training loop, having the solver and the automatic differentiation stack in the same language and the same process is the point, and that is a stronger reason than any individual algorithm.

What the readme does not give you is a benchmark table. There is a benchmark directory in the repository, and the readme links to the DiffEqTutorials.jl repository for the Henon-Heiles example, but the choice between adaptive methods for a specific problem is left to the documentation and to your own measurement, which is the right way round for a package with this many solvers.

Editorial conclusion

Adopt OrdinaryDiffEq.jl if you are solving differential equations in Julia and want the full solver menu of the SciML ecosystem behind one import, particularly if you need symplectic or split methods for Hamiltonian systems or event handling for discontinuities. Do not adopt it casually on an existing project, because v7 carries breaking changes across all sublibraries, and the readme directs you to DifferentialEquations.jl if you want the batteries-included interface. Before you upgrade, read NEWS.md for the migration guide, and if you are coming from Python, weigh the cost of moving your whole numerical stack rather than just the solver.

Frequently asked questions

Official sources

1. Issues
2. Project website
3. README
4. SciML/OrdinaryDiffEq.jl on GitHub