% This ODE stepper takes a single step of length h, breaking it into as many
% sub-steps as necessary to achieve the desired absolute or fractional
% tolerance. The step doubling is implemented with recursion.
% CCK 2009-Feb-25
function [x1, greatest_rec_depth] = adaptive(ODE, x0, t0, h, stepper, order,... 
   tolerance, maximum_recursion_depth, recursion_depth)

if ~exist('tolerance', 'var')
   tolerance = [1e-5, 1e-5];
end

if ~exist('recursion_depth', 'var') 
   recursion_depth = 1;
end

if ~exist('maximum_recursion_depth', 'var') 
   maximum_recursion_depth = 8;
end


% tolerance has the form [abs, rel]
abs_tol = tolerance(1);
rel_tol = tolerance(2);

if (recursion_depth > maximum_recursion_depth) 
   abs_tol
   rel_tol
   maximum_recursion_depth
   recursion_depth
   time = t0
   error('Maximum recursion depth exceeded:')
end

% first call to the step doubler
[x1, err] = twostep(ODE, x0, t0, h, stepper, order);

greatest_rec_depth = recursion_depth;
   % This will change if the step doubles again below!
   
% If *both* relative and absolute tolerances are exceeded, then 
% do it in two steps.
if (  ( abs(max(err)) > abs_tol ) &&...
      ( norm(err)/norm(x1-x0) > rel_tol )  )

   [xhalf, NR1] = adaptive(ODE, x0, t0, h/2, stepper, order, tolerance/2,...
      maximum_recursion_depth, recursion_depth + 1 );
      
   [x1, NR2] = adaptive(ODE, xhalf, t0+(h/2), h/2, stepper, order, tolerance/2,...
      maximum_recursion_depth, recursion_depth + 1 );
   % Note the recursive calls to adaptive!!!
   
   greatest_rec_depth = max(NR1, NR2);   
end

