% This is a front end for ODE-solving using adaptive.m. It is a slight
% improvement of adaptive_stepper --- allows the user to specify all
% the desired time points explicitly (like Octave's lsode).
% CCK 2009-Mar-01
function [x, greatest_depth] = adaptive_stepper2(ode, x0, t,... 
   stepper, order, tolerance, max_rec_depth)

Nt = numel(t); % Number of time points, including initial conditions
Ndim = numel(x0); % How big is the phase space?
x = zeros(Nt, Ndim); % each row is a state at a given timestep.
x(1,:) = x0;

greatest_depth = 0;  % initialize greatest depth
for n = 2:Nt
   h = t(n) - t(n-1);
   [x(n,:), depth] = adaptive(ode, x(n-1,:), t(n-1), h, stepper, order,...
      tolerance, max_rec_depth); 
      % N.B.: pass function handle as ode rather than @ode.
   greatest_depth = max(depth, greatest_depth);
end