Let me get to the saving part last.

The documentation for the callbacks is here:

https://juliadiffeq.github.io/DiffEqDocs.jl/latest/features/callback_functions.html

It’s a little underdocumented because it’s going to be changed so that way it’s easier to make it to every solver on the common interface (currently it’s only implemented for OrdinaryDiffEq.jl. It’s not hard to add it to others, what’s hard is making it somewhat uniform while keeping the features).

If you’re just modifying the array values of the solution, you can just do:

```
my_callback = @ode_callback begin
for i = 1:length(f.L)
u[i] *= ...
end
@ode_savevalues
end
```

The changing of all of the cache is only needed if you’re changing non-index values or if you’re changing the size of the ODE system on the fly (maybe there are other odd cases like that).

However, note that if you are changing `u`

like this at every step, the solution may no longer be continuous. This means that anything relying on interpolations (`saveat`

and `dense`

) will no longer work properly. If you do need this, you may want to have this just be an event which occurs every step, which is just by making `event_f(args...) = true`

and then making an `apply_event!`

. The downside here is that, to appropriately handle the discontinuity it has to save both before and after the event application.

May I ask why you are doing this post modification? If it’s to set some kind of normalization (multiplication, so that looks like what you may be doing?) or other condition on `u`

, you probably should be using a DAE solver instead.

Pre-allocating in Julia doesn’t do too much because `push!`

is amortized. The reason why I can make this statement and know it for a fact is that this part of the interface exists, it’s just undocumented. If you actually look at the definition for `solve`

in OrdinaryDiffEq.jl, you will see that it is:

```
function solve{uType,tType,isinplace,T<:OrdinaryDiffEqAlgorithm,F}(
prob::AbstractODEProblem{uType,tType,isinplace,F},
alg::T,timeseries=[],ts=[],ks=[];<then it defines the kwargs>
```

so you can pass in `timeseries`

, `ts`

, and `ks`

where `timeseries`

is a `Vector`

of the solutions at each timepoint, `ts`

are the timepoints, and `ks`

is a vector for saving interpolation information (no need here since `dense=false`

). Note that it has to match what will be put in there, so `timeseries`

is a `Vector`

where each value matches the type of `u0`

. So if `u0`

is a vector of size 2^18, `timeseries`

should be a `Vector{typeof(u0)}(num_timepoints)`

. Note that this is the same setup as the solution, so you can just check `sol.u`

on a similar problem if you aren’t sure.

I actually spent quite some time making this “safe”, as in it doesn’t even need to be the length of the final solution. It will fill the pre-allocated vector until it’s full, and then `push!`

into it to finish (or if it doesn’t make it that far, it will `resize!`

at the end to match the true length). So you don’t actually have to worry about making it the right size, just make sure it’s

But note that this isn’t implemented in other solvers yet. Actually, report back on whether it makes a difference. If it makes a difference on your system, maybe that means I just didn’t test it hard enough. Whether or not it makes a difference, please let me know. If it does make a substantial difference, open up an issue and we’ll get this implemented in the other solvers as well.

Now lastly:

The setup you want to use is the following:

```
solve(prob, Tsit5(), abstol=1e-10, dense=false,
save_timeseries=false, tstops=zs,
saveat = zs)
```

Maybe the `save_timeseries`

argument should be renamed. Your setup (`save_timeseries=true`

and `tstops`

used) means is save every internal step of the method. `tstops`

means that the solvers will stop at each value in `tstops`

. So what this means is that it will try to step as large as possible (getting the tolerance), except if in that interval it has a value of tstops, in which case it will shorten the interval to hit `tstops`

exactly, and save every step.

What the setup I am giving you means: "step as far as possible so that you hit the correct tolerance or a value in `tstops`

, and ‘interpolate’ to save the values only at the timepoints in `zs`

". Note from earlier interpolation won’t actually work because your solution is discontinuous at each solution step, but since we make `tstops`

the same as all of the values in `saveat`

, there is no interpolation that is needed.

That said, I think you actually hit a bug: the interpolation won’t actually save anything if hits the value directly. I put a fix up on OrdinaryDiffEq.jl master and will put a patch out for this tonight.

There has been a discussion to change the argument to `save_timeseries`

to `save_everystep`

or something like that because this has come up in the chat before. If you want to pursue something like this, please open an issue. Also, since this saving behavior is modular (mix and matching `tstops`

, `saveat`

, `save_timeseries`

, `dense`

, and `calck`

), I should probably make a tutorial notebook showing all of the different combinations and how to make it do things like this. (That would be a great contribution to DiffEqTutorials.jl if anyone finds the time!).