SavingCallback influences result of ODE problem

Specific Domains Modelling & Simulations
10

Hi Chris, sorry for the delay. The code I provided indeed doesn’t run, but is just an indicator of where the problem occurred. Here is a working example:

using DifferentialEquations, DiffEqCallbacks, SparseArrays, LinearAlgebra

function build(Nn::Int64,BC::String)::Matrix{Float64}
    if BC=="open"
        Jm=spdiagm(1=>ones(Nn),-1=>ones(Nn)); Jm[1,2]*=2;
    elseif BC=="periodic"
        Jm=spdiagm(1=>ones(Nn),-1=>ones(Nn)); Jm[1,2]*=2;
        Jm[Nn+1,Nn+1]=1.0;
    elseif BC=="neumann"
        Jm=spdiagm(1=>ones(Nn),-1=>ones(Nn)); Jm[1,2]*=2;
        Jm[Nn+1,Nn]*=2;
    else
        error("No matching boundary condition")
    end
    return Jm
end
function initialChain(n0::Int64,Nn::Int64,ϕ::Float64,IC::String)::Vector{ComplexF64}
    ψ0=zeros(ComplexF64,Nn+1)
    if IC=="full" # Central site is full
        ψ0[:].=sqrt(n0)*exp(im*ϕ);
    elseif IC=="empty" # Central site starts of empty
        ψ0[2:Nn+1].=sqrt(n0)*exp(im*ϕ);
        ψ0[1]=sqrt(0.1*n0)*exp(im*ϕ);
    else
        error("Undefined initial condition")
    end
    return ψ0
end
function du!(du,u,p,t)                  # Differential equations
    J,Jm,U,γ,D=p
    mul!(D,Jm,u)
    @. du = im*J*D - im*U*abs2(u)*u
    du[1]+= -0.5γ*u[1]
end

function affect_ic!(integrator)         # Callback functions
    integrator.p[1]+=integrator.p[end]
end
save_func(u,t,integrator)=integrator.p[1]
saved_values=SavedValues(Float64,Float64)

U=1; n0=700; Nn=200; J=70; γ=150; ϕ=0.0pi; IC="full"; BC="periodic";
tObs=Array(0:0.001:1); dosetimes=Array(range(0.5/U/40,0.5/U,40))
D=zeros(ComplexF64,Nn+1); Jm=build(Nn,BC); ψ0=initialChain(n0,Nn,ϕ,IC)

if IC=="full"
    δJ=(J-100.0)/40
    param=[100.0, Jm, U, γ, D, δJ]
elseif IC=="empty"
    δJ=(J-40.0)/40
    param=[40.0, Jm, U, γ, zeros(ComplexF64,M), δJ]
end
cb_preset=PresetTimeCallback(dosetimes,affect_ic!,save_positions=(false,false))
cb_save=SavingCallback(save_func,saved_values,saveat=tObs)
cb=CallbackSet(cb_preset,cb_save)

prob=ODEProblem(du!, ψ0, (tObs[1], tObs[end]), param, callback=cb)
sol=solve(prob,Vern8(),reltol=1e-8,abstol=1e-8,adaptive=true,saveat=tObs)


using Plots
plot(tObs,abs2.(sol[1,:])/n0)

Using only the PresetTimeCallback, that changes p, works fine. Taking lazy=false does not seem to change the result. Using the SavingCallback does require lazy=false, even when you leave out the PresetTimeCallback.