Code breaks on an oscillator-like problem

Hello,

I have a system of two coupled ODEs. I tried several solvers from Diffrax, both non-stiff and stiff, with and without time step control. Every solver falls on this problem with the same error:

equinox.EquinoxRuntimeError: The maximum number of solver steps was reached. Try increasing `max_steps`.

However, they do solve the system when setting max_steps=1_000_000. The funny thing is that this is a test calculation to find a solution in a small vicinity of a center point of the ODE system, i.e. I expected the small amplitude oscillations around the initial values.

I suspect that I'm doing something wrong with Diffrax, because the scopy.integrate.solve_ivp does my system with DOP853 method. So the jax.experimental.ode.odeint does it too.

Here, the snippet:

import jax.numpy as jnp
import diffrax

def rhs(t, w, a):
    x, y = w
    xdot = - x**2 + y**2 - a**2
    ydot = - 2*x*y
    return xdot, ydot

a = 16.
w0 = (0., 1.1*a)

tbeg, tend, dt = 0., 1e3, 1.
saveat = diffrax.SaveAt(ts=jnp.arange(tbeg, tend, dt))

term = diffrax.ODETerm(rhs)
solver = diffrax.Dopri8()
stepsize_controller = diffrax.PIDController(rtol=1e-8, atol=1e-8)
#stepsize_controller = diffrax.PIDController(rtol=1e-3, atol=1e-6,
#                                            pcoeff=0.4, icoeff=0.3, dcoeff=0)
sol = diffrax.diffeqsolve(term, solver, tbeg, tend, 1e-3*dt, w0, args=a,
                          saveat=saveat, stepsize_controller=stepsize_controller)

Start by running the inline Diffrax reproducer with the reported solvers, tolerances, and max_steps behavior. Compare the solver step count and resulting trajectory against the stated SciPy and JAX baselines, then trace the relevant Diffrax step-size or termination behavior. Done means identifying whether the failure is a Diffrax bug or an expected configuration issue and documenting a reproducible fix or explanation.