I want to track all variables within Julia’s DifferentialEquation’s efficiently (easy and fast).
My current solution, on a test function, is:

function test2(du,u,p,t)
global i
global Out_values
global sampletimes
a,b,c = p
d=a^0.1*(t+1)
e=u[1]/a
f=u[2]/d
if t in sampletimes
Out_values[1,i] = d
Out_values[2,i] = e
Out_values[3,i] = f
i=i+1
end
du[1] = a*u[1]
du[2] = d*u[2]
du[3] = b*u[2] - c*u[3]
end
sampletimes = tspan[1]:0.3:tspan[2]
Out_values = Array{Float64}(undef, 3, 2*11)
i=1
prob = ODEProblem(test2,u0,tspan,p)
sol = solve(prob,saveat=0.3,tstops=sampletimes)

However, there are two problems:

Somehow I get duplicated results (not sure why) with slightly different results, and

I can’t track variables if I decide to use all results, i.e. sol = solve(prob).

A second approach could also be to compute the additional values afterwards. For example you could define d_values = @. p.a^(0.1)*(sol.t), e_values = @. sol[1,:] ./ p.a and f_values = @. sol[2,:] ./ p.d… If you don’t have too many values to track, this might be easier.

To answer the questions:

You get duplicated results because solvers (especially implicit solvers) evaluate the function many times, to approximate the Jacobian, or for adaptive stepsizes. These additional evaluations are tracked in your approach. The callback library takes care of these things!

I don’t understand the second problem fully, but with the saving callbacks you could call the solve function with any arguments (i.e. also without the saveat). With the second approach it directly uses the sol array, so it shouldbe compatible in all cases.

@DrPapa, thank you very much for answering my question. I took a look at the package, and it seems it serves the purpose. I will investigate more later and check how it performs against the callback library that @SteffenPL proposed.

@SteffenPL, thanks for the answer. Is there a tutorial or an example of a similar problem where callbacks are applied?

It’s for example in the documentation when you scroll down.

For you case it looks like this:

function saving_callback(u, t, integrator)
a, b, c = integrator.p
d=a^0.1*(t+1)
e=u[1]/a
f=u[2]/d
return (d, e, f)
end
saved_values = SavedValues(Float64, Float64, Float64)
cb = SavingCallback(saving_callback, saved_values)
sol = solve(prob, callback=cb)

now saved_values.saveval contains the values corresponding to the times saved_values.t.