Hello

I have some code for simulating complex differential equations of a certain model that used to work fine, but fails after going from Julia 1.7 to 1.8 and updating all the packages.

```
using DifferentialEquations, DiffEqCallbacks, SparseArrays
function du!(du,u,p,t) # Define differential equations
J,Jm,U,γv,D=p
mul!(D,Jm,u)
@. du = im*J*D - (im*U*abs2(u)+γv)*u
end
save_func(u,t,integrator)=integrator.p[1]
saved_values=SavedValues(Float64,Float64)
cb_save=SavingCallback(save_func,saved_values,saveat=tObs)
M=size(ψ0,1); mat=build(jp); # Set parameters
param=[jp.J, mat.Jm, jp.U, mat.γv, zeros(ComplexF64,M)]
#---------------------------------------
prob=ODEProblem(du!, ψ0, (tObs[1], tObs[end]), param, callback=cb_save) # Define ODE problem & solve
sol=solve(prob,Vern8(),dt=dt0,adaptive=true,saveat=tObs)
```

Without the SavingCallback there is no problem and the `Vern8()`

solver is more than suffucient with the default tolerances. When the Savingcallback is added, the exact same setup gives incorrect results. In the figure I plot the absolute value of one element of the vector `u`

that in this model is supposed to stabilize.

I find it counter-intuitive that simply storing a parameter changes the solving algorithm. Am I missing something trivial here?