function [tout, yout] = rk45(FunFcn, t0, tfinal, y0, tolrel,tolabs)
%RK45	Integrate a system of ordinary differential equations using
%	the classical 4th order Runge-Kutta formula.
%
%	[T,Y] = RK45('yprime', T0, Tfinal, Y0) 
%       integrates the system of ordinary differential equations 
%	described by the M-file YPRIME.M over the interval T0 to Tfinal 
%	and using initial conditions Y0.
%	[T, Y] = RK45(F, T0, Tfinal, Y0, TOLREL,TOLABS) uses relative and
%       tolerances TOLREL and TOLABS. The relative tolerance is measured with 
%	respect to the norm of the initial value Y0
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution column-vector.
%         yprime - Returned derivative column-vector; yprime(i) = dy(i)/dt.
% t0    - Initial value of t.
% tfinal- Final value of t.
% y0    - Initial value column-vector.
% tolrel- The desired relative accuracy. (Default: tolrel = 1.e-3).
% tolabs- The desired absolute accuracy. (Default: tolabs = 1.e-3).
%
% OUTPUT:
% T  - Returned integration time points (column-vector).
% Y  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, yout).


% Initialization

% default tolerance
% at t=tfinal the norm of the error should be smaller than tolabs 
% thus for a step of size h the error should be of the order 
%       tolabs * h/(tfinal-t0)
if nargin < 5, tolrel = 0.001; end
if nargin < 6, tolabs = 0.001; end
tol = max([tolabs,tolrel*norm(y0)])/(tfinal-t0);

% at least 5 steps, at most 10000 steps
t = t0;
hmax = (tfinal - t)/5;
hmin = (tfinal - t)/10000;
h = (tfinal - t)/100;
y = y0(:);
tout = t;
yout = y.';


% The main loop
while  (t < tfinal) 
  if t + h > tfinal, h = tfinal - t; end
    
  % Compute the slopes
  % one Runge-Kutta step of size h
    k1 = feval(FunFcn, t, y);
    k2 = feval(FunFcn, t+h/2, y+h*k1/2);
    k3 = feval(FunFcn, t+h/2, y+h*k2/2);
    k4 = feval(FunFcn, t+h,   y+h*k3);
    res1 = y + (k1 + 2*k2 + 2*k3 + k4)*h/6;
    
  % two Runge-Kutta steps of size h/2
    h=h/2;
    % use k1 from above 
    k2 = feval(FunFcn, t+h/2, y+h*k1/2);
    k3 = feval(FunFcn, t+h/2, y+h*k2/2);
    k4 = feval(FunFcn, t+h,   y+h*k3);
    res2 = y + (k1 + 2*k2 + 2*k3 + k4)*h/6;
    k1 = feval(FunFcn, t+h,     res2);
    k2 = feval(FunFcn, t+3*h/2, res2+h*k1/2);
    k3 = feval(FunFcn, t+3*h/2, res2+h*k2/2);
    k4 = feval(FunFcn, t+2*h,   res2+h*k3);
    res2 = res2 + (k1 + 2*k2 + 2*k3 + k4)*h/6;
    h = h*2;

  % Estimate the error
    yerr=norm(res1 - res2);

  % stepsize too large
    if ( yerr > (h*tol) && (h > hmin))
      h = 0.8*h;
    else
     % accept values, do extrapolation
     % use the LOCAL error estimates
      t = t+h;
      y = (16*res2-res1)/15;
      tout = [tout; t];
      yout = [yout; y'];
      % make h larger if possible
      if ( yerr< 0.5*h*tol)
        h = min(hmax, 1.2*h);
      end
    end
  end;

  if (t < tfinal)
    disp('SINGULARITY LIKELY.')
    t
  end
