function [x,y,stats] = ode4545(dfun,xspan,y0,AbsTol,RelTol,thresh)
%*************************************************************************
% [x,y,stats] = ode4545(dfun,xspan,y0,c,AbsTol,RelTol)
%
% Integrates ODE system from xspan(1) to xspan(2) by the
% explicit (4,5)-pair Runge-Kutta (RK) formula with stage order 3 and (in
% TeX notation) "$\hat{b}_7 \not= 0$",
% automatically regulating step size so that
%
%    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   Error tolerances
%     [optional]
%   thresh [optional] Threshold for changing step size between steps
%   
% Outputs:
%   x               Mesh of algorithmically-chosen points in xspan
%   y               Solution of DE (5th-order approximation) on the mesh x
%   stats           Optional output 
%                   [# steps, # failed steps, # function evals]
% 
% Demo (requires lnhde.m):
%   xspan=[0, 2*pi]; y0=0;
%   [x,y] = ode4545('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 Runge-Kutta')
%   figure(2), plot(x,yex-y)
%   title(sprintf('Maximum error is %8.2e',norm(yex-y,'inf')))
%
% Demo (requires orbit.m):
%   xspan = [0, 6.19216933]; y0 = [1.2; 0; 0; -1.04935750983031990726];
%   [x,y] = ode4545('orbit',xspan,y0);
%   plot(y(1), y(2));
%
% Demo (requires arenstorf.m):
%   xspan = [0 17.0652165601579625588917206249];
%   y0 = [0.994; 0; 0; -2.00158510637908252240537862224];
%
% Demo (requires kepler.m):
%   Set ecc to some number between 0 and 1 (eccentricity of orbit)
%   xspan = [0 2*pi]; y0 = [1 - ecc; 0; 0; ((1+ecc)/(1-ecc))^(1/2)]
%*************************************************************************

if nargin < 6           % Provide default values for optional parameters
    thresh = 0.5;
  if nargin < 5
    RelTol = 1e-5;
    if nargin < 4
      AbsTol = 1e-4;
    end
  end
end

%---------------------------------------------------------------

% Coefficients
c1 = 0; c2 = 1/5; c3 = 3/10; c4 = 4/5; c5 = 8/9; c6 = 1; c7 = 1;
bh1 = 71/57600; bh2 = 0; bh3 = -71/16695; bh4 = 71/1920;
bh5 = -17253/339200; bh6 = 22/525; bh7 = -1/40;
a31 = 3/40; a32 = 9/40; a41 = 44/45; a42 = -56/15; a43 = 32/9;
a51 = 19372/6561; a52 = -25360/2187; a53 = 64448/6561; a54 = -212/729;
a61 = 9017/3168; a62 = -355/33; a63 = 46732/5247; a64 = 49/176;
a65 = -5103/18656;
b1 = 35/384; b2 = 0; b3 = 500/1113; b4 = 125/192; b5 = -2187/6784;
b6 = 11/84;

%---------------------

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 number of dependent variables
x = [x0; zeros(32,1)];    % Create mesh
y = [y0'; zeros(32,neq)]; % Allocate space for solution
init = 1;                 % Flag: starting integration

% Stats
numfailed = 0;
numderiv = 0;

% MAIN LOOP: Take steps one at a time until end of interval reached.

while dr*(xf-x) > 0
  if init

% First Step: Find first initial stepsize using
% simple estimates of |df/dx| and |df/dy|

    dy0 = feval(dfun,x0,y0);

    delta = sqrt(eps)*max([abs(x0) 1]);
    normdfx = norm((feval(dfun,x0+delta,y0)-dy0) / delta, inf);
    deltay = sqrt(eps)*max([abs(y0') 1]);
    normdfy = norm((feval(dfun,x0,y0+deltay)-dy0) / ...
              norm(deltay, inf), inf);
    normf = norm(dy0, inf);
    
    h = sqrt(2)*eps^(1/6) / sqrt(normdfx + normdfy*normf);
    h = min([h abs(xf-x0)/8]); % h is at most 1/8 of integration interval
    h = dr*h;
    hfn = h*dy0;
    
    numderiv = numderiv + 3;
    
    init = 0;
    
  else

%   Not first step.
%   Select step based on error estimate in step just completed.

    up = (thresh*(AbsTol+RelTol*norm(y(n,:),inf)) / Est)^(1/6);
    h = up*h;
    hfn = up*hfn;
  end

% Attempt fifth-order Runge-Kutta step and estimate error
 
  k1    = hfn;
  k2    = h*feval(dfun,x(n)+c2*h,y(n,:)'+ c2*k1);
  k3    = h*feval(dfun,x(n)+c3*h,y(n,:)'+(a31*k1+a32*k2));
  k4    = h*feval(dfun,x(n)+c4*h,y(n,:)'+(a41*k1+a42*k2+a43*k3));
  k5    = h*feval(dfun,x(n)+c5*h,y(n,:)'+(a51*k1+a52*k2+a53*k3+a54*k4));
  k6    = h*feval(dfun,x(n)+c6*h,y(n,:)'+(a61*k1+a62*k2+a63*k3+a64*k4 ...
          +a65*k5));
  k7    = h*feval(dfun,x(n)+c7*h,y(n,:)'+(b1*k1+b2*k2+b3*k3+b4*k4 ...
          +b5*k5+b6*k6));
  ynp1  = y(n,:)' + (b1*k1+b2*k2+b3*k3+b4*k4+b5*k5+b6*k6);  % 1st step

  numderiv = numderiv + 6;

  Est = norm(h*(k1*(b1-bh1) + k3*(b3-bh3) + k4*(b4-bh4) + k5*(b5-bh5) ...
        + k6*(b6-bh6) + k7*(0-bh7)),inf);
      
  if Est < AbsTol+RelTol*norm(y(n,:),inf)

%   Passed error test: accept step, advance solution

    x(n+1) = x(n) + h;
    y(n+1,:) = ynp1';
    n = n+1;
    hfn = k7;
    
    if mod(n,32)==0                % Allocate more space for solution
      x = [x; zeros(32, 1)];
      y = [y; zeros(32, neq)];
    end
  else
      
    numfailed = numfailed + 1;
    
  end
end

x1 = [x(n-2), x(n-1), x(n)];
y1 = [y(n-2, :)', y(n-1, :)', y(n, :)'];
dy = [feval(dfun,x(n-2),y(n-2,:)'), feval(dfun,x(n-1),y(n-1,:)'),...
     feval(dfun,x(n),y(n,:)')];

numderiv = numderiv + 3;
    
yco = HermInt(x1, y1, dy);

[m1,n1] = size(yco);
for k = 1:m1
    sol(k) = Horner2(x1,yco(k,:),xf);
end
    
% Record number of computational steps prior to truncation
[numsteps one] = size(x);

% Assigns sol and truncates x depending on xf
  
if xf <= x(n-1)
  x(n-1) = xf;
  y(n-1, :) = sol;
  x = x(1:n-1);   % Chop off unused array space
  y = y(1:n-1, :);
else
  x(n) = xf;
  y(n, :) = sol;
  x = x(1:n);     % Chop off unused array space
  y = y(1:n, :);
end

% Either outputs stats as third parameter or as text
if nargout > 2
    stats = [numsteps;numfailed;numderiv];
else
    % Display stats
    disp(sprintf('%d successful steps', numsteps))
    disp(sprintf('%d failed attempts', numfailed))
    disp(sprintf('%d function evaluations', numderiv))
end

end