function [tout, yout] = rk4(FunFcn, tspan, y0, ssize)

% RK4	Integrates a system of ordinary differential equations using
%	the fourth order Runge-Kutta  method.  See also ODE45 and
%	ODEDEMO.M.
%	[T,Y] = RK4('yprime', Tspan, Y0) integrates the system
%	of ordinary differential equations described by the M-file
%	YPRIME.M over the interval Tspan=[t0,tfinal] and using initial
%	conditions Y0.
%	[T, Y] = RK4(F, Tspan, Y0, SSIZE) uses step size SSIZE
%
% 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.
% tspan = [t0, tfinal], where t0 is the initial value of t, and tfinal is
%         the final value of t.
% y0    - Initial value column-vector.
% ssize - The step size to be used. (Default: ssize = (tfinal - t0)/100).
%
% 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(T,Y).


% Initialization

t0=tspan(1);
tfinal=tspan(2);

if nargin < 4, ssize = (tfinal - t0)/100; end

h=ssize;

% Initialization
t = t0;
y = y0(:);
dt = abs(tfinal - t0);
N = floor(dt/ssize) + 1;
if (N-1)*ssize < dt
  N = N + 1;
end
tout = zeros(N,1);
tout(1) = t;
yout = zeros(N,size(y,1));
yout(1,:) = y.';
k = 1;

% The main loop
while (k < N)
  if t + h > tfinal 
    h = tfinal - t; 
    tout(k+1) = tfinal;
  else
    tout(k+1) = t0 +k*h;
  end
  k = k + 1;

  % Compute the slopes 
  x=y; 
   thold=t; yhold=y; 
  s1 = eval(FunFcn); s1 = s1(:);
  t=thold+h/2; y=yhold + h*s1/2;
  x=y;
  s2 = eval(FunFcn); s2=s2(:);
  t=thold+h/2; y=yhold + h*s2/2;
  x=y;
  s3 = eval(FunFcn); s3=s3(:); 
  t=thold+h; y=yhold + h*s3;
  x=y;
  s4 = eval(FunFcn); s4=s4(:);
   t=thold; y=yhold; 
  y = y + h*(s1 + 2*s2 + 2*s3 +s4)/6;
  t = tout(k);
  yout(k,:) = y.'; 
  if max(yout(k,:))==inf 
	disp(['Step size ',num2str(h),' gives infinite result for t = ',num2str(t)]) 
	break 
  end 		
end;