Different results for dopri5 depending on loop vs solve

I am trying to solve the Lorenz butterfly using diffrax, and in particular the Dopri5 solver.

I take two approaches.

In of them, I simply solve over the entire time range (0,40), using a step size DISCRETE_DT=1e-2

In another, I solve over each time range (t, t+DISCRETE_DT), and then loop this solver.

Annoyingly, the numerical error between these 2 approaches grows with the sequence length?

Does anyone have a good explanation for this phenomenon? I know the Lorenz butterfly is chaotic and that Dopri5 takes adaptive step sizes but

• I'm using float64
• I've checked and no steps are being rejected in either trajectory. So, I would have thought the exact same numerical computations are being used in both settings.

Here is a colab with a repro: https://colab.research.google.com/drive/1aNjUCB4SmQ_UDJ3H5Lqu4jzfFxOwY6Q6?usp=sharing

Start with the linked Colab reproduction and compare the full-range solve with the looped per-interval Dopri5 solves using float64. Trace how the solver chooses or carries step sizes across these two approaches, and determine why their numerical results diverge. Done means the cause is explained clearly and any needed follow-up is identified.