I'm interested in having a deval function. Since odeXX is native to Octave now and does have the output form of sol = odeXX(...), it seems like what's there should be enough to extract deval from the existing functions: It already interpolates internally.

The ode solvers in Matlab can be called with two main different
forms for the output

or

the second form does not merely aggregate outputs in one
structure but triggers a completely different behaviour of the
solver.

The idea is that with the second form the solver advances until
the end of the timespan (or passed it) and does not care of
evaluation of the state at intermediate timesteps.

The output structure "out" contains all data that can be fed
to "deval" to evaluate the solution at any time point between
tspan(1) and tspan(end) to the maximum possible order of
accuracy consistent with the solution method used.

We don't have deval in Octave at the moment and I cannot find
documentation of the actual fields of the output structure in
matlab, so I think we don't need to copy those fields exactly.

On the other hand here is the list of data we would need to do
interpolation as a post-processing step (at least for explicit runge kutta):

times selected by the adaptive selection algorithm and state evaluated there