I can’t try this right now, but I don’t understand why you have save_positions = (false,false). Have you run it without that?

I’ve never used that option - this is what the manual says:

save_positions : Boolean tuple for whether to save before and after the affect! . The first save will always occcur (if true), and the second will only occur when an event is detected. For discontinuous changes like a modification to u to be handled correctly (without error), one should set save_positions=(true,true) ."

Makes no difference for the result. However, if set to (true, true) I get double entries in the time vector of the solution.
I also don’t get why a correct computation would depend on this setting.

I don’t either but read that line in the docs and thought it was worth a try. I’m going to try to run your example once some things finish installing and precompiling.

ODE integrators are usually making a polynomial approximation to the true solution. If you have a change in slope at the event, there is no polynomial approximation for that. The (true, true) tells the integrator to use different approximations for each side of the event. It no longer feels compelled to make things C1 continuous.

By the way, I would also advise against a discrete event such as if 10. <= int.t <= 20. 1., which gives the integrator no continuous information to triangulate the event. Better to provide a zero-crossing like the bouncing ball demo in docs. A continuous function helps the integrator know when to adjust its steps as it gets closer to the event.

Also by the way, if you do a discrete event, you needn’t do if ... 1 else 0, you can just make it logical, 10. <= int.t <= 20.. But again, better to do a continuous zero crossing.

Thank you for the input but I do not think that this is the issue here. Otherwise why would the version of @klaf using the parameter in the callback work?
In the bouncing ball example the problem is that the event is not detected. When using the periodic callback you already know the time of the event so there should be no problem.

When using the DEDataArray, you have to remember to update all of the cache variables. Notice the docs show the full_cache iterator.

When solving an ODE with discontinuities, you have to make sure you resolve the discontinuities if you want a high precision answer. Here, you know them in advance, so you can set tstops = [10.0,20.0] to hit those points exactly during the integration.

With those two facts together it all works out:

using OrdinaryDiffEq, Plots
mutable struct SimType{T} <: DEDataVector{T}
x::Vector{T}
a::T
end
f1(du, u, p, t) = du[1] = u.a
f2(du, u, p, t) = du[1] = if 10. <= t < 20. 1. else 0. end
function cb!(int) # the callback function
for c in full_cache(int)
c.a = if 10.0 <= int.t < 20. 1. else 0. end
end
return
end
u0 = SimType([0.], 0.)
tspan = (0., 30.)
prob1 = ODEProblem(f1, u0, tspan)
sol1 = solve(prob1, Tsit5(); callback = PeriodicCallback(cb!, 1., save_positions = (false,false)))
prob2 = ODEProblem(f2, [0.], tspan)
sol2 = solve(prob2, Tsit5(), tstops = [10.,20.])
function trueSol(t)
if t <= 10.
0.
elseif 10. < t <= 20.
-10. + t
else
10.
end
end
plot(sol1.t, trueSol.(sol1.t), lw=2, label="true solution")
plot!(sol1.t, [x[1] for x in sol1.u], linestyle = :dash, legend = :topleft, lw=2, label="solution with callback")
plot!(sol2.t, [x[1] for x in sol2.u], linestyle = :dot, lw=5, label="solution without callback")

It definitely can’t be used with Sundials, and generally with implicit it will have undefined behavior. I’m probably removing it from the docs because of these issues, and I wouldn’t quite recommend it anymore.