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

% usage:  [tout, yout] = rk4(FunFcn, t0, tfinal, 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', T0, Tfinal, Y0) integrates the system
%	of ordinary differential equations described by the M-file
%	YPRIME.M over the interval T0 to Tf and using initial
%	conditions Y0.
%	[T, Y] = RK4(F, T0, Tfinal, 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.
% t0    - Initial value of t.
% tfinal- 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(tout, yout).

% Written by John C. Polking, Rice University
% Last modifed February 3, 1995

% Initialization

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

h=ssize;

% Initialization
t = t0;
y = y0(:);
tout = t;
yout = y.';

% The main loop
while (t < tfinal)
      	if t + h > tfinal, h = tfinal - t; end

% Compute the slopes
  	s1 = feval(FunFcn, t, y); s1 = s1(:);
      	s2 = feval(FunFcn, t + h/2, y + h*s1/2); s2=s2(:);
	s3 = feval(FunFcn, t + h/2, y + h*s2/2); s3=s3(:);
	s4 = feval(FunFcn, t + h, y + h*s3); s4=s4(:);
	t = t + h;
	y = y + h*(s1 + 3*s2 + 2*s3 +s4)/6;
	tout = [tout; t];
	yout = [yout; y.'];

end;

