%  Example of the use of Rke and Yvalue in MATLAB.
%  REQUIRES Rke.m, Yvalue.m, and X2rkef.m.

%  Solves van der Pol's equation on the interval [0,20] and plots the 
%  answers obtained at each step as well as several values obtained by 
%  interpolation in each step.  In a phase plane plot, the extra values 
%  are needed for a smooth graph.  

%  Set initial values and tolerances.
x = 0;
lastx = 20;
y = [1; 1];
tol = 1e-5;
threshold = [1e-5; 1e-5];

%  Initialize output arrays to initial values.    
X = 0;
Y = y';

h = 0.1;
nasteps = 0;
yp = zeros(2,1);
while x <= lastx
  [x,y,yp,h,nasteps,flag,step,ycoeff] = ...
    Rke('X2rkef',x,y,yp,h,tol,threshold,nasteps);
  if flag ~= 0
    error(['flag = ',int2str(flag),' at x = ',num2str(x),'.']);
  end

  %  Compute several approximations in [x-step,x] to get a smooth graph.
  %  Do not produce approximations past lastx:
  for i = 3:-1:0
    xi = min(x,lastx) - i*(step/4);
    z = Yvalue(xi,x,step,ycoeff);
    X = [X; xi];
    Y = [Y; z'];
  end
end

plot(X,Y)   
legend('y(x)','y''(x)')
title('van der Pol equation')

figure
plot(Y(:,1),Y(:,2))
title('phase plane plot')
xlabel('y(x)')
ylabel('y''(x)')
   
---Cut here Rke.m
function [x,y,yp,h,nasteps,flag,step,ycoeff] = ...
                       Rke(f,x,y,yp,h,tol,threshold,nasteps)
%  Rke integrates a system of neq first order ordinary differential 
%  equations, y' = f(x,y), over one step using a Runge-Kutta method
%  due to R. England. Yvalue can be used after any successful step by
%  Rke to approximate the solution components in the interval [x-step,x].
%
%  Input arguments:
%   f     = name of the function to evaluate y' = f(x,y).  The function
%           should have the form
%              function 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.
%   yp    = column vector yp = f(x, y).  On the first call its value does
%           not matter, but on subsequent calls, it is to be the value
%           returned by Rke.
%   h     = step size for current step.  The sign of h 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.  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 for each i
%              abs(local error in y(i)) <= tol*max(threshold(i),abs(y(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 returned by Rke.
%
%  Output arguments:
%   x      = the integration has advanced to this value of the independent
%            variable.
%   y      = column vector of computed solution values at x.
%   yp     = column vector yp = f(x,y).
%   h      = step size suitable for the next step.
%  nasteps = number of steps attempted.
%   flag   = 0  for a successful step.
%            1  if tol or threshold is too small.
%
%  Optional output arguments:
%   step   = output x - input x, the actual step taken.
%   ycoeff = array of coefficient values used by Yvalue.
%
%  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 by setting nasteps = 0.
%      flag = 1  the error tolerances 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 tolerances, set
%                nasteps = 0 and call Rke again.

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

%  Test some input data.
if (tol > 0.01) | (tol < 10 * eps)
  error('tol must satisfy 10*eps <= tol <= 0.01.')
end

tolaim = 0.6*tol;
first = false;

if nasteps == 0

  %  First call: test more input data.
  first = true;
  valid = true;
  if any(threshold < 0)
    valid = false;
  elseif any((abs(threshold)+abs(y)) <= 0)
    valid = false;
  end
  if ~valid
    error('threshold has an illegal value.')
  end
       
  %  Initialize yp.
  yp = feval(f, x, y);
end

failed = false;
while true    %  Loop until an h is found for which the step succeeds.

  %  Compute a 4th order result ymid at xmid = x + h/2.
  k1 = feval(f, x+0.25*h, y + 0.25*h*yp);
  k2 = feval(f, x+0.25*h, y + 0.125*h*(yp + k1));
  k3 = feval(f, x+0.5*h, y + 0.5*h*(-k1 + 2*k2));
  xmid = x + 0.5*h;
  ymid = y + (h/12)*(yp + 4*k2 + k3);
  k4 = feval(f, xmid, ymid);

  %  Compute a 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 a 5th order result y5th at x + h.
  ytemp = y + (h/12) * (-yp - 96 * k1 + 92 * k2 - 121 * k3 + ...
              144 * k4 + 6 * k5 - 12 * k6);
  ktemp = feval(f, x+h, ytemp);
  y5th = y + (h/180) * (14 * yp + 64 * k2 + 32 * k3 - 8 * k4 + ...
              64 * k6 + 15 * k7 - ktemp);

  %  Form err, the weighted maximum norm of the estimated local error.
  wt = (2 * (abs(y) + abs(ymid)) + abs(y4th) + abs(y5th)) / 6;
  wt = max(wt, threshold);
  delta = abs(y5th - y4th);  %  delta(i) is an estimate of the error in y4th(i).
  positive = find(wt > 0);   %  wt(i) = 0 => delta(i) = 0.
  delta(positive) = delta(positive) ./ wt(positive);
  err = norm(delta,inf);

  nasteps = nasteps + 1;
  if err > tol

    %  Failed step.
    if first
      h = 0.1*h;
    elseif ~failed
      failed = true;
      h = max(0.1,(tolaim/err)^0.2)*h;
    else
      h = 0.5*h;
    end
    if abs(h) < 24*eps*max(abs(x), abs(x+h))
      flag = 2;
      return
    end
  else

    %  Successful step.
    x = x + h;
    k8 = feval(f, x, y5th);
    if nargout >= 7
      step = h;
    end
    if nargout == 8
      alpha = h * (yp - k8);
      beta = y5th - 2 * ymid + y;
      gamma = y5th - y;
      ycoeff(:,1) = ymid;
      ycoeff(:,2) = h * k4;
      ycoeff(:,3) = 4 * beta + 0.5 * alpha;
      ycoeff(:,4) = 10 * gamma - h * (yp + 8 * k4 + k8);
      ycoeff(:,5) = -8 * beta - 2 * alpha;
      ycoeff(:,6) = -24 * gamma + 4 * h *(yp + 4 * k4 + k8);
    end
    y = y5th;
    yp = k8;
    temp = 1 / max(0.1, (err/tolaim)^0.2);
    if failed
      temp = min(1, temp);
    end
    h = temp * h;
    flag = 0;
    return
  end

end