Let's talk details then. On the small stiff ODE and non-stiff end SUNDIALS is bested quite easily (for example https://benchmarks.sciml.ai/html/StiffODE/Hires.html), but that's pretty well known so not much to say there. For the larger ODE end, it's missing Runge-Kutta-Chebyschev methods, exponential integrators, etc. which do amazing on a lot of equations, and in direct head to head FBDF has been outperforming CVODE BDF for over a year now on most equations (as demonstrated in the online benchmarks, and now there are some applications papers starting to show that). Where SUNDIALS does do well is for distributed NVector is really good, along with it's heterogeneous compute. Its new multirate methods are good too. But then it's lacking methods for SDEs, DDEs, RODEs, higher precision, probabistic numerics, etc. so it's hard to call it feature filled. It's also lacking a complete symbolic front end like ModelingToolkit for handling higher index DAEs, so while IDA currently outperforms DFBDF head-to-head, the optimal thing to do is still to use a Julia solver just post MTK optimized (for example see https://benchmarks.sciml.ai/html/DAE/ChemicalAkzoNobel.html and for the large equation side see some of the stuff on HVAC models). Sundials has deep customizability in the linear and nonlinear solvers but that's also seen on the Julia side.

And while SUNDIALS does have an adjoint implementation, it does not construct the vjp through reverse mode so direct (naive) usage is really lacking in performance. Yes I know that there is a way to define the vjp function in the docs, but it's buried enough in there that all of the examples I've looked at of SUNDIALS in applications in the wild (like in AMIGO2) don't make use of it so there's a real performance loss seen by end users for many cases. DiffEq of course automates this.

Don't get me wrong, SUNDIALS is really good, still bests the Julia solvers on the distributed case, and there is a small performance advantage to CVODE in the one doc example on the 2D Brusselator after adding iLU preconditioning (https://diffeq.sciml.ai/dev/tutorials/advanced_ode_example/ and the low tolerance case of the Filament benchmark (https://benchmarks.sciml.ai/html/MOLPDE/Filament.html, though the medium tolerance case is bested by ROCK4). For years DifferentialEquations.jl used wrappers to SUNDIALS as a crutch for the cases where SUNDIALS still outperformed it, so nowadays it both has a fairly complete wrapper of SUNDIALS but also has benchmarks of different algorithms outperforming it.

We've been benchmarking against it weekly for years now, and within the last two years we've really seen a reversal so now it's a rare sight to see SUNDIALS perform the best. And there's a lot of momentum too: new parallel solvers, new GPU solvers, etc. are coming out from both crews so there's a fun arm's race going on right now, so we will see how it evolves (I definitely want to wrap those multirate methods ASAP and play with the theory a bit)

> Btw. several "research" papers at machine learning conferences on "Neural Differential Equations" presented "discoveries" covered in the documentation of Sundials https://sundials.readthedocs.io/en/latest/cvodes/Mathematics....

Of course the Julia folks know this because we were the ones mentioning this at conferences, since we would routinely choose to benchmark against SUNDIALS implementations instead of random unstable ML ones.

Independently on our views on the comparisons of Sundials and DifferentialEquations.jl, I want to wholeheartedly agree with you that modern research in Neural Differential Equations and ML in general is rediscovering a ton of stuff that is obvious of everyone that has done serious numerical work.