% This is a front end for ODE-solving using adaptive.m.
% CCK 2009-Feb-25
function [x,t, greatest_depth] = adaptive_stepper(ode, x0, t0, t1, Nsteps,... 
   stepper, order, tolerance, max_rec_depth)


h = (t1 - t0)/Nsteps;
t = (t0 + [0:Nsteps]*h)'; % column vector of times
Ndim = numel(x0); % How big is the phase space?
x = zeros(Nsteps, Ndim); % each row is a state at a given timestep.
x(1,:) = x0;

greatest_depth = 0;  % initialize greatest depth
for n = 2:Nsteps+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