function [rkfT, rkfW, rkfH] = rkfsys(yprime, t0, tlast, w0, hmax, hmin, TOL)
%function [rkfT, rkfW, rkfH] = rkfsys(yprime, t0, tlast, w0, hmax, hmin, TOL)
%%%%%%% RUNGE - KUTTA - FEHLBERG  ALGORITM FOR SYSTEMS  %%%%%%%%%%%%%%%%%%
%% 
%%  Requires problem set up in an m-file yprime: w0 is a vector of
%%  initial conditions, output rkfW is a matrix, each column
%%  corresponding to  a solution at one of the equations
% INPUT:
% yrime   - String containing name of problem m-file (enclose in '  ')
% t0    - Initial value of t.
% tlast - Last value of t.
% w0    - Vector of initial values.
% hmax  - maximum step size;  defaults to (tlast-t0)/10
% hmin  - minimum step size;  defaults to (tlast-t0)/200
% TOL   - local truncation tolerence; defaults to 1.e-6
%
% OUTPUT:
% rkfT  - time points (column-vector).
% rkfW  - Solution column-vectors per rkfT-values.
% rkfH  - Column of step sizes h
% The result can be displayed by: plot(rkfT, rkfW).

if nargin < 4, 
 error('Not enough arguments, enter help rkf53_2 for help')
end;
if nargin < 5, 
 hmax=(tlast-t0)/10;
end;
if nargin < 7, 
 TOL = 1.e-6;
end;
if nargin < 6, 
 hmin=(tlast-t0)/200;
end;
% Initialization
rkfT = t0;
h = hmax;
rkfW = w0;
rkfH=0;
flag=1;   
t=t0;w=w0;
     %% Step 2   the main loop
while flag == 1;
     %% Step 3 the function evaluations
 k1=h*feval(yprime,t,w);
 k2=h*feval(yprime,t+h/4,w+k1/4);
 k3=h*feval(yprime,t+(3*h/8),w+(3*k1/32)+(9*k2/32));
 k4=h*feval(yprime,t+(12*h/13),w+(1932*k1/2197)-(7200*k2/2197)+(7296*k3/2197));
 k5=h*feval(yprime,t+h,w+(439*k1/216)-8*k2+(3680*k3/513)-(845*k4/4104));
 k6=h*feval(yprime,t+(h/2),w-(8*k1/27)+2*k2-(3544*k3/2565)+(1859*k4/4104)-...
      (11*k5/40));
     %% Step 4 the error control function R=|wtil-w|/h;
  R=max(abs((k1/360)-(128*k3/4275)-(2197*k4/75240)+(k5/50)+(2*k6/55)))/h;
     %% Step 5   computing the step size factor
  delta=0.84*(TOL/R)^(1/4);
     %% Step 6  Check tolerence
    if R<=TOL    
        t=t+h;   %% if OK, compute next approximation  Step 7
        w=w+(25*k1/216)+(1408*k3/2565)+(2197*k4/4104)-k5/5;
      rkfT=[rkfT;t];rkfW=[rkfW;w]; rkfH=[rkfH;h];  %% Step 8  OUTPUT result
    end;
      %% Step 9   Adjust step size
    if delta <=0.1;
         h=0.1*h;   %% not too small
    elseif delta >= 4;
         h=4*h;    %% not too big
    else
         h=delta*h;  %% just right
    end;
       %% Step 10   check h against hmax
    if h > hmax;
       h=hmax;
    end;
        %% Step 11  Check stopping criteria
    if t >= tlast;
        flag=0;
    elseif   t+h > tlast;
        h=tlast-t;
    elseif    h < hmin;
        flag=0;
        message='minimum step size violated'
    end;
end;
