function [x,y,k0,h,flag,step,ycoeff,nasteps] = ...
                       Rke(f,x,y,k0,h,tol,threshold,nasteps)
%  Rke integrates a system of neq first order ordinary differential 
%  equations over one step using a Runge-Kutta method due to R. England.
%  It provides for the evaluation of the solution within the step by
%  interpolation.  The quintic Hermite interpolant is evaluated in Yvalue.  
%  Yvalue can be used independently after any successful step by Rke to 
%  approximate the solution components in the interval [x-step,x].
%
%  Input parameters:
%   f     = a string giving the name of the function defining the system 
%           of differential equations.  The function must have the form
%                yprime = f(x, y)
%           where yprime is a column vector.
%   x     = initial value of the independent variable.
%   y     = column vector of neq solution values at x.
%   k0    = column vector k0 = f(x, y).  On the first call its value does
%           not matter, but it must be present as a column vector of the
%           same length as y.  On subsequent calls, it is to be the value
%           output by Rke.
%   h     = step size for current step (its sign determines the
%           direction of integration).  On the first call to Rke you
%           must specify h.  After the first call, the code suggests
%           a suitable h for the next step.
%   tol, threshold = desired tolerances: tol is a scalar and threshold
%           is a column vector with neq components.  The value of tol must 
%           be in the interval [ 10*eps , 0.01 ] where eps is the unit
%           roundoff for your machine; each component of threshold must
%           be non-negative.  On the first call, if some y(i) = 0, then
%           threshold(i) must be positive for that i.  The convergence
%           criterion is roughly
%                abs(local error in y(i)) <=
%                      tol*max(threshold(i),abs(y(i)))
%           for each i.
%   nasteps = number of steps attempted.  On the first call to Rke you
%           must set nasteps = 0.  On subsequent calls, it is to be the
%           the value output by Rke.
%
%  Output parameters:
%   x      = value of the independent variable to which the
%            integration advanced.
%   y      = column vector of computed solution values at x.
%   k0     = column vector k0 = f(x, y).
%   h      = step size suitable for the next step.
%   flag   = an integer reporting what the code did:
%              0, if a step was successfully taken;
%              1, if excessive work (more than 500 steps) was expended;
%              2, if tol, threshold are too small;
%   step   = output x - input x, the actual length of step taken.
%   ycoeff = array of coefficient values for the quintic Hermite 
%            interpolation used by Yvalue.
%   nasteps = number of steps attempted.
%
%  Rke is organized so that subsequent calls to continue the integration 
%  involve little, if any, extra effort.  The value of flag must, however, 
%  be monitored in order to determine what to do next.  Specifically, if
%      flag = 0, the code may be called again to continue the
%                integration another step towards x+h.  You may
%                shorten h if you wish; if its sign is changed,
%                you must restart the integration (nasteps = 0).
%      flag = 1, the code has attempted 500 steps: either tol/threshold 
%                are not appropriate or possibly the solution does not 
%                exist beyond x.  You may alter tol or threshold.  If you
%                want to continue, set nasteps = 0 and call Rke again.
%      flag = 2, the error parameters tol/threshold are too small.
%                If a solution component is zero, use a positive value
%                for the corresponding component of threshold.  To
%                continue with increased error parameters, set
%                nasteps = 0 and call Rke again.

false = 0; true = ~false;
maxsteps = 500;
neq = length(y);
ycoeff = zeros(neq,6);

if nasteps == 0

  %  First pass; test the input data.
  valid = (tol <= 0.01) & (tol >= 10 * eps);
  negative = find(threshold < 0);
  valid = valid & isempty(negative);
  yzero = find(y == 0);
  if ~isempty(yzero)
    both_zero = find(threshold(yzero) == 0);
    valid = valid & isempty(both_zero);
  end
  if ~valid
    error('tol or threshold has a value not allowed in Rke.')
  end
       
  %  Initialize.
  k0 = feval(f, x, y);
  first = true;
else

%  Continuation; test some input data.
  first = false;
  if (tol > 0.01) | (tol < 10 * eps)
    error('tol has a value not allowed in Rke.')
  end
end

%  Try to take a step.    
tolaim = 0.6*tol;
failed = false;
flag = -2;
while flag == -2
  %  Compute 4th order result ymid at xmid = x+h/2.
  k1 = feval(f, x+0.25*h, y + 0.25*h*k0);
  k2 = feval(f, x+0.25*h, y + 0.125*h*(k0 + k1));
  k3 = feval(f, x+0.5*h, y + 0.5*h*(-k1 + 2*k2));
  ymid = y + (h/12)*(k0 + 4*k2 + k3);
  xmid = x + 0.5*h;
  k4 = feval(f, xmid, ymid);

  %  Compute 4th order result y4th at x+h.
  k5 = feval(f, xmid+0.25*h, ymid + 0.25*h*k4);
  k6 = feval(f, xmid+0.25*h, ymid + 0.125*h*(k4 + k5));
  k7 = feval(f, xmid+0.5*h, ymid + 0.5*h*(-k5 + 2*k6));
  y4th = ymid + (h/12)*(k4 + 4*k6 + k7);

  %  Compute 5th order result y5th at x+h.
  tempy = y + (h/12) * (-k0 - 96 * k1 + 92 * k2 - 121 * k3 + ...
              144 * k4 + 6 * k5 - 12 * k6);
  tempk = feval(f, x+h, tempy);
  y5th = y + (h/180) * (14 * k0 + 64 * k2 + 32 * k3 - 8 * k4 + ...
              64 * k6 + 15 * k7 - tempk);

  %  Form err, the weighted maximum norm of the estimated local error.
  wt = (2 * (abs(y) + abs(ymid)) + abs(y4th) + abs(y5th)) / 6;
  wt = max(threshold, wt);
  delta = abs(y5th - y4th);
  positive = find(delta > 0);
  if isempty(positive)
    err = 0;
  else
    err = norm(delta(positive) ./ wt(positive),inf);
  end
  if err > tol
    %  Failed step.
    if first
      alpha = 0.1;
    elseif ~failed
      failed = true;
      alpha = (tolaim/err)^0.2;
      alpha = max(0.1, alpha);
    else
      alpha = 0.5;
    end
    h = alpha * h;
    if abs(h) < 24*eps*max(abs(x), abs(x+h))
      flag = 2;
    end
  else
    %  Successful step.
    x = x + h;
    k8 = feval(f, x, y5th);
    step = h;
    alpha = step * (k0 - k8);
    beta = y5th - 2 * ymid + y;
    gamma = y5th - y;
    ycoeff(:,1) = ymid;
    ycoeff(:,2) = step * k4;
    ycoeff(:,3) = 4 * beta + 0.5 * alpha;
    ycoeff(:,4) = 10 * gamma - step * (k0 + 8 * k4 + k8);
    ycoeff(:,5) = -8 * beta - 2 * alpha;
    ycoeff(:,6) = -24 * gamma + 4 * step *(k0 + 4 * k4 + k8);
    y = y5th;
    k0 = k8;
    first = false;
    beta = (err/tolaim)^0.2;
    alpha = 1 / max(0.1, beta);
    if failed
      alpha = min(1, alpha);
    end
    h = alpha * h;
    flag = 0;
  end
  nasteps = nasteps + 1;
  if (nasteps == maxsteps)
    flag = 1;
  end 
end

