Hey everyone!!

I am stuck on trying to replicate some previous results for my thesis right now.

For that, I want to generate the solution to a system of differential equations using the OrdinaryDiffEq package.

Previously, those results were generated using Julia 1.5.3. Now I am using Julia 1.10.

For some reason I am unable to replicate them even when using the exact same lines of code.

My system of differential equation is a kuramoto model, solved using Tsit5 (Runge-Kutta). Afaik (not a physicist though), this solver and model should be numerically stable and deterministic, aka yield the same results for the same set of initial conditions and parameters.

Here are the relevant lines of code:

```
"""
kuramoto_model(ddu, du, u, kura_param, t)
The Kuramoto model of coupled oscillators.
# Arguments
- `ddu`: The second derivative of the state vector.
- `du`: The first derivative of the state vector.
- `u`: The state vector.
- `kura_param`: The parameters of the Kuramoto model.
- `t`: The current time.
# Notes
This function modifies `ddu` in-place.
"""
function kuramoto_model(ddu, du, u, kura_param, t)
V = exp.(u .* im)
flow = real.(V .* conj(kura_param.L * V))
@. ddu = -kura_param.kuram_Î± .* du + kura_param.P + flow
end
kura_problem = SecondOrderODEProblem(kuramoto_model, du0, u0, tspan, kura_param_data) # initial ODE problem
sol = solve(kura_problem, Tsit5(); save_idxs = 1:(d.N), reltol = 1e-6, abstol = 1e-6, callback = cb, saveat = save_at_range)
```

Neither callback nor save_at_range changes between the julia versions.

When I evaluate this code using 1.5.3 I get a slightly different result than when I evaluate it using 1.10.

These changes can not be attributed to floating point errors I think, as the maximum of the state variable at the last timestep sometimes differs by more than 1.

Examples:

julia 1.5.3 yields `max(sol.u[end]) = 4.166`

while 1.10 yields `max(sol.u[end])=3.3866`

The initial conditions and all kuramoto parameters are the same for both versions, i triple-checked that.

Does anyone here know if something changed in between those julia versions that could explain this drastic difference? Am I overlooking something simple?

Thank you very much in advance, Iâ€™ve been stuck on this bug (?) or user error for the past two weeks and donâ€™t know how to continue. If you need further information lmk, Iâ€™ve created an entire powerpoint trying to debug this lol (not kidding).

hereâ€™s an example plot of one of the solutions where 1.5 yields 0.0822 and 1.10 yields 0.0031: