diffeqpy: DifferentialEquations.jl solvers called from Python

A Python front end over the Julia solver stack

diffeqpy is a package for solving differential equations in Python, and it utilizes DifferentialEquations.jl for its core routines to give high performance solving of many different types of differential equations. The list of equation types is the reason to choose it. It handles discrete equations including function maps and discrete stochastic simulations such as Gillespie and Markov processes, ordinary differential equations, split and partitioned ODEs including symplectic integrators and IMEX methods, stochastic ordinary differential equations, random differential equations, differential algebraic equations, delay differential equations, and mixed discrete and continuous equations such as hybrid equations and jump diffusions. The mechanism underneath is Julia running inside the Python process, and the setup.py makes that explicit, since the install requirements are juliacall at 0.9.28 or newer and jill. pyjulia is named in the README as the bridge with its own documented limitations, and juliacall is the package that provides it. So this is not a reimplementation of these solvers in Python, it is a way to reach the Julia implementations from Python, and that is the whole design.

de. or ode., and why the choice matters

Setup is a single import, and there are two of them depending on how much you need.

from diffeqpy import de

The de name brings in the solvers available in DifferentialEquations.jl. If only the solvers available in OrdinaryDiffEq.jl are required, the README gives the alternative.

from diffeqpy import ode

The difference is practical rather than semantic, since ode. has a shorter loading time and a smaller memory footprint compared to de. Most commands work without any modification, and the general flow is to follow exactly what would be done in Julia except to add de. or ode. in front. That design has a direct consequence for how you learn the tool, because the DifferentialEquations.jl documentation and the DiffEqTutorials are the main in-depth documentation for this package, and the README spends its length showing how to translate those Julia examples into Python code. If you already know the Julia ecosystem, this package feels familiar immediately. If you do not, the learning curve is the Julia solver API rather than a Python-shaped one, and the README is upfront that the upstream docs are where the detail lives.

The bang problem and the _b convention

One wrinkle in the translation is worth understanding before you hit it, because it is a naming collision with Python itself. The Julia API uses a trailing bang on mutating functions, so step! mutates the integrator in place, and Python does not allow a ! in function names, which the README notes is a limitation of pyjulia. The convention is that functions which have a ! on the Julia side are written with _b instead.

from diffeqpy import de

def f(u,p,t):
    return -u

u0 = 0.5
tspan = (0., 1.)
prob = de.ODEProblem(f, u0, tspan)
integrator = de.init(prob, de.Tsit5())
de.step_b(integrator)

This is the integrator interface rather than the solve interface, and it is the mode you use when you want to step a solution manually, for example to inspect state between steps or to implement your own control loop. The de.init call takes the problem and an algorithm, here the Tsit5 method, and returns the integrator object that step_b then advances. The renaming is not cosmetic. Any example you copy from Julia documentation will have bangs in it, and a search for the function will fail in Python unless you translate the name, so the _b suffix is the thing to internalise before you start.

Solving an ODE and reading the result

The basic solve follows the Julia shape, with a Python derivative function, an initial condition, a time span, and a solve call.

from diffeqpy import de

def f(u,p,t):
    return -u

u0 = 0.5
tspan = (0., 1.)
prob = de.ODEProblem(f, u0, tspan)
sol = de.solve(prob)

The solution object is the same as the one described in the DiffEq tutorials and in the solution handling documentation, with one caveat the README flags in parentheses, that the array interface is missing. The practical fields are sol.t for the solution time points and sol.u for the solution values, and the interpolation sol(t) gives a continuous solution, which is how you get a finer grid without re-solving.

import numpy
t = numpy.linspace(0,1,100)
u = sol(t)
plt.plot(t,u)
plt.show()

The missing array interface is the one thing to check against your own code, because code that expects sol.u to behave like a numpy array will not get that, and the plotting example in the README has to stack the values explicitly for a system of equations. The common interface arguments control the solve, so a method choice and output interval and tolerances are all passed to the solve call.

sol = de.solve(prob,de.Vern9(),saveat=0.1,abstol=1e-10,reltol=1e-10)

de.jit and seval, and the limits of compilation

Because the derivative function sits in the solver's inner loop, its speed matters, and the README offers two ways to improve it. The first is de.jit, which converts the entire ODE problem to symbolic form, optimizes that symbolic form and emits efficient native code to simulate it, improving the efficiency of the solver at the expense of added setup time.

fast_prob = de.jit(prob)
sol = de.solve(fast_prob)

The second option is to define the function in Julia directly, which also allows for specialization and could be helpful for repeat or long calls.

jul_f = de.seval("(u,p,t)->-u") # Define the anonymous function in Julia
prob = de.ODEProblem(jul_f, u0, tspan)
sol = de.solve(prob)

The limitation is stated plainly. de.jit uses ModelingToolkit.jl's modelingtoolkitize internally and some restrictions apply, and not all models can be jitted, with the README pointing at the modelingtoolkitize documentation for the details. That is the catch worth planning around, because the symbolic path only accepts models that can be expressed in the form the tool recognises, and a model with unusual structure, non-standard functions or awkward control flow may simply not be jittable. The seval route has no such restriction, since it hands the function to Julia directly, but it moves your model out of Python.

Systems of ODEs, and the stack and transpose helpers

For a system of equations the pattern is the same with an array initial condition and an array-valued derivative, and the Lorenz system is the example the README works through.

def f(u,p,t):
    x, y, z = u
    sigma, rho, beta = p
    return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]

u0 = [1.0,0.0,0.0]
tspan = (0., 100.)
p = [10.0,28.0,8/3]
prob = de.ODEProblem(f, u0, tspan, p)
sol = de.solve(prob,saveat=0.01)

plt.plot(sol.t,de.transpose(de.stack(sol.u)))
plt.show()

The parameters are passed through the problem object as a third argument to the constructor, matching the Julia signature where the derivative takes parameters, and the initial condition is a plain Python list. The two helper calls in the plotting line are the part to notice. de.stack turns the list of solution vectors into something array-shaped, and de.transpose handles the orientation, since the solver stores solutions as a sequence of state vectors and the plotting call wants them the other way round. The same stack call is what you would use for a phase plot, where you index the three components of the stacked array directly.

us = de.stack(sol.u)
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(us[0,:],us[1,:],us[2,:])
plt.show()

The 3D plot also needs the Axes3D import even in recent matplotlib versions, which is a small compatibility wart in the example rather than in the library.

The Julia runtime, the licence and the release record

The dependency question dominates every practical decision about this package, because install_requires is juliacall and jill, which means installing diffeqpy brings a Julia runtime into the process. jill is the installer that fetches Julia, so a first install pulls a working Julia toolchain, and that download and the startup cost of initialising the Julia session are the real cost of the approach. The package metadata describes the environment as Python 3.10 and later, with classifiers up to 3.14, and python_requires set to 3.10 or newer, so the Python side is not the constraint. The licence is MIT, and the authors are Chris Rackauckas and Takafumi Arakaki, with the contact address at the SciML organisation and the keywords listing the equation classes along with oneapi, AMD, ROCm, CUDA and Metal, which points at the GPU backends available on the Julia side. The release record shows v2.6.0 on 2026-05-14, v2.5.5 on 2026-01-21 and v2.5.0 in December 2024, with the last push on 2026-10-17 and the repository not archived. Two Colab notebook examples are linked, one on solving the Lorenz equation faster than SciPy and Numba and one on solving ODEs on GPUs, which are the fastest way to see whether the approach suits your problem before you commit to installing a Julia toolchain locally.

Editorial conclusion

Adopt diffeqpy if you already work in Python but need a solver choice or an equation class that SciPy does not cover, and you are willing to have a Julia runtime in the process. Skip it if you need a pure Python install with no Julia toolchain, or if a SciPy solver is enough, since the extra machinery buys you the Julia ecosystem rather than speed for simple problems. Verify first that your platform can install the Julia runtime that juliacall and jill pull in, that your derivative function can be expressed in a form the JIT path accepts if you plan to use de.jit, and that the array interface limitation on the solution object does not break your plotting or analysis code.

Frequently asked questions

Official sources

1. Issues
2. License: MIT
3. README
4. Releases
5. SciML/diffeqpy on GitHub