function [x,y] = EulerAuto(dfun,xspan,y0,AbsTol,RelTol)
% [x,y] = EulerAuto(dfun,xspan,y0,AbsTol,RelTol)
%
% Integrates ODE system from xspan(1) to xspan(2) by
% Euler method, automatically regulating step size so
%
%    ErrEst < AbsTol + RelTol*norm(y,inf).
%
% Inputs:
%   dfun    string that names a function M-file for derivatives
%   xspan   two-component vector -- initial & final values of 
%           independent variable
%   y0      initial value(s) of dependent variable(s)
%   AbsTol,RelTol [optional] error tolerances
%   
% Outputs:
%   x       mesh of automatically chosen pts in xspan
%   y       solution of DE (Euler approximation) on the mesh x
% 
% Demo (requires lnhde.m):
%   xspan = [0 2*pi]; y0=0;
%   [x,y] = EulerAuto('lnhde',xspan,y0);
%   yex   = (25/13)*cos(x) + (5/13)*sin(x) - (25/13)*exp(-5*x);
%   figure(1), plot(x,y,'b',x,yex,'r')
%   title('Solution of dy/dx = -5y + 10 cos(x) by Euler Formula')
%   figure(2), plot(x,yex-y)
%   title(sprintf('Maximum error is %8.2e',norm(yex-y,inf)))

if nargin < 5
  RelTol = 1e-3;
  if nargin < 4
    AbsTol = 1e-6;
  end
end

x0 = xspan(1);             % Initial point
xf = xspan(2);             % Final point
dr = sign(xf-x0);          % Direction of integration
n  = 1;                    % Number of steps completed
y0 = y0(:);                % Column vector
[neq one] = size(y0);      % Determine # of dep. vars.
x  = [x0; zeros(32,1)];    % Create mesh
y  = [y0'; zeros(32,neq)]; % Allocate space for solution
f0 = feval(dfun,x0,y0);    % Initial derivative
nde   = 1;                 % Count derivative evals
nfail = 0;                 % Count failed steps
init  = 1;                 % boolean: starting integration
 
% MAIN LOOP: Take steps one at a time until end of interval reached.
 
while dr*(xf-x) > 0
  if init
 
%   First Step: Estimate 2nd deriv bound with
%    steps 1/8, 1/64, 1/512 interval
 
    halfM2 = sqrt(eps);
    for p = 1:3
      h      = (xf-x0)/8^p;
      f1     = feval(dfun,x0+h,y0+h*f0);
      nde    = nde + 1;
      E      = norm( (h/2)*(f1-f0), inf);
      halfM2 = max([halfM2 E/h^2]);
    end
 
%   Select initial stepsize -- choose a step that
%   is expected to pass the error test, but is at
%   most 1/8 of the integration interval.
   
    Tol = AbsTol+RelTol*norm(y0,inf);
    h = sqrt(0.8*Tol/halfM2);
    h = min([h abs(xf-x0)/8]);
    h = dr*h;
    hfn = h*f0;
    init = 0;
  else
 
%   Not first step.
%   Select step based on error estimate in step just completed.
   
    Tol = AbsTol+RelTol*norm(y(n,:),inf);
    up  = sqrt(0.8*Tol/Est);
    h   = up*h;
    hfn = up*hfn;
  end
 
% Attempt an Euler step, estimate error
 
  ynp  = y(n,:)' + hfn; 
  hfnp = h*feval(dfun,x(n)+h,ynp);
  nde  = nde+1;
  Est  = 0.5*norm(hfnp-hfn,inf);
 
  if Est < AbsTol+RelTol*norm(y(n,:),inf)
 
%   Passed error test -- accept step, advance solution
%   [use local extrapolation and evaluate derivative
%    at the corrected point]
   
    ynp    = y(n,:)' + 0.5*(hfn+hfnp);
    n      = n+1;
    x(n)   = x(n-1) + h;
    y(n,:) = ynp';
    hfn    = h*feval(dfun,x(n),ynp);
    nde    = nde+1;
    if mod(n,32)==0     % Allocate more space for solution
      x = [x; zeros(32,1)];
      y = [y; zeros(32,neq)];
    end
  else
    nfail = nfail+1;
  end
end
x = x(1:n);   % Integration complete; delete unused array space
y = y(1:n,:);

% Interpolate to deliver solution at final point



% Print Stats
disp(sprintf('%d successful steps',n))
disp(sprintf('%d failed attempts',nfail))
disp(sprintf('%d function evaluations',nde))