Maybe you can find some useful solvers in this repo:
The solvers are compatible with solve_ivp, in case you need to use event handling.
You can also go “halfway” using the Python wrappers of DifferentialEquations.jl: GitHub - SciML/diffeqpy: Solving differential equations in Python using DifferentialEquations.jl and the SciML Scientific Machine Learning organization · GitHub
It does come with the limitation,
Autodiff does not work on Python functions. …
so you’d need to move (at least) your callback functions to pure Julia if autodiff is a must.
There is rklib (GitHub - jacobwilliams/rklib: Fixed and variable-step Runge-Kutta solvers in Modern Fortran · GitHub) which has a big collection of variable-step explicit Runge-Kutta methods.
Otherwise I’m not aware of anything comparable in terms of composability of different ODE solvers, linear solvers, non-linear solvers, event handling facilities and other bells and whistles. The investment in the SciML libraries is enormous. On the other hand we have LLMs nowadays so I believe it’s feasible to create your own set of SciPy-compatible or custom set ODE solvers, as long as you’ve already got your autodiff library figured out (PyTorch, JAX, …?).
I have investigated a number of ODE solving methods in the past (see Time-stepping method for cylindrical diffusion equation). Maybe if you provide more information about your problem we can give a better recommendation.
Here are few points to consider:
- size of the ODE system: 10, 100, 1000, …
- stiff or non-stiff?
- differential-algebraic equations
- mass matrix
- event handling
- symbolic Jacobian, numerical differencing, autodiff…
- dense, banded, sparse, or other Jacobian pattern
- discrete PDE problem
I’ve written Python bindings to the RKC ODE solver in Fortran before (https://github.com/ivan-pi/rkc/tree/python/python). It took me ~1 day of work and this was without LLMs. I did have a decent understanding of C/Python/Fortran interfacing, which made things easier. What I found the most challenging was learning to use Meson and figuring out how to create a pip-compatible Python Package.
When it comes to speed, the scipy.integrate.OdeSolver expects the solver to perform a single-step. This is so the abstract parent class can drive the integration and do the event-handling (if used). But it is not optimal from a performance point of view, because it means at each ODE solver step you go the full route Python → C → Fortran → C → Python, unpacking and re-packing the argument tuples, switching from interpreted to compiled code and vice-versa, which adds overhead.
To avoid this overhead you’d want to write your own solve_ivp like driver, and do all the time-stepping in the compiled language, ideally using a callback function that is also compiled (either calling C/Fortran or using Numba). I noticed this issue when designing the RKC bindings (https://github.com/ivan-pi/rkc/blob/python/python/rkc/rkc_proc.py), but I didn’t get around to “pushing” the time loop down into the C layer yet.