This is a follow-up question on DiffEq skips over event in the first `solve` but works after `reinit!` - #16 by ChrisRackauckas. The code simulates a circuit in which an AC power source is supplying an RL load through a diode. I used a `SavingCallback`

to save `u[1]`

at fixed intervals, but the result contains data point(s) that are extremely large.

The code below is a bit long, but the suspect issue is related to the `SavingCallback`

.

```
using DifferentialEquations
using Plots
using BenchmarkTools
function I_load!(du, u, p, t)
du[1] = (sqrt(2)*230*sin(2*pi*50.0*t) - u[1] * p[1]) *2
end
dt = [];
# test if diode needs to be turned off
function condition_off(u, t, integrator)
# println("test affect off at $(integrator.t), u=$(integrator.u)")
# println("In condition_off")
# @show integrator.t integrator.dt
# push!(dt, integrator.dt)
if integrator.p[2] == 1
u[1]
else
return 1
end
end
function affect_off!(integrator)
integrator.p[1] = 1e10 # total resistance (load + diode)
integrator.p[2] = 0
# println("In affect_off!")
end
# test if diode needs to be turned on
function condition_on(u, t, integrator)
# println("In condition_on")
if integrator.p[2] == 0
# voltage across the diode minus the minimum turn-on voltage (1)
ret = (sqrt(2)*230*sin(2*pi*50.0*t) - (0.5 * get_du(integrator)[1] + u[1] * 100)) - 1
else
ret = -1
end
end
function affect_on!(integrator)
# println("applied affect on at t=$(integrator.t)")
integrator.p[1] = 100 # total resistance (load + diode)
integrator.p[2] = 1
end
cb_off = ContinuousCallback(condition_off, nothing, affect_off!,
repeat_nudge=1//10) # only when condition goes negative
cb_on = ContinuousCallback(condition_on, affect_on!, nothing,
repeat_nudge=1//10) # only when condition goes positive
saved_values = SavedValues(Float64, Tuple{Float64})
function save_u(u, t, integrator)
@show integrator.t u[1] # <- outputs come from here
(u[1], )
end
cb_save = SavingCallback(save_u, saved_values,
saveat=0:1/2000:0.06,
)
cb_set = CallbackSet(cb_off, cb_on, cb_save);
u0 = [0.0]
p = [100., 1.] # resistance, diode flag
tspan = (0.0, 0.06)
prob = ODEProblem(I_load!, u0, tspan, p)
solver = Rodas5P()
integrator = init(prob, solver,
callback=cb_set,
abstol=1e-5,
reltol = 1e-6,
dtmax=1/60,
saveat=0:1/2000:0.06;
);
sol = solve!(integrator)
plot(sol, label="I Rodas5P")
```

Output:

```
integrator.t = 9.999999999999999e-5
u[1] = 0.0
integrator.t = 0.0010355149302289026
u[1] = 0.02466443779613447
integrator.t = 0.0010355149302289026
u[1] = 0.09489608425988974
integrator.t = 0.0016705724725223623
u[1] = 0.20457147552261043
integrator.t = 0.002348954342345277
u[1] = 0.3470141746597498
integrator.t = 0.0030822014110460615
u[1] = 0.5151192104287117
integrator.t = 0.0030822014110460615
...truncated...
integrator.t = 0.052876422347891805
u[1] = 0.49957763366843355
integrator.t = 0.053380394153775476
u[1] = -1.0304261248426776e40 # <- this does not seem correct
integrator.t = 0.05444431272273977
u[1] = -3.528390922450684e-8
integrator.t = 0.05444431272273977
u[1] = -2.0560975872248562e-8
integrator.t = 0.06
u[1] = -3.213201469371499e-8
integrator.t = 0.06
u[1] = -3.254155700933845e-8
integrator.t = 0.06
u[1] = -3.2109302361073413e-8
integrator.t = 0.06
u[1] = -3.0880836404674164e-8
integrator.t = 0.06
u[1] = -2.890823031148805e-8
integrator.t = 0.06
u[1] = -2.625004076845552e-8
integrator.t = 0.06
u[1] = -2.2971309978110092e-8
integrator.t = 0.06
u[1] = -1.9143565658578353e-8
integrator.t = 0.06
u[1] = -1.4844821043579996e-8
integrator.t = 0.06
u[1] = -1.0159574882427828e-8
integrator.t = 0.06
u[1] = -5.1788114400275675e-9
integrator.t = 0.06
u[1] = -4.968781781758675e-16
```

The issue is that it silently fails. If the `abstol`

and `reltol`

are set to other values, the extremely large values will appear at other times (but all close to discontinuities).