% function y = rk4(a,b,n,y0,f)
%  
%  Runge-Kutta method for systems
%
%    y'  = f(t,y),  t in [a,b]
%    y(a) = y0
%
%  Inputs
%   a,b    time interval
%   n      number of subintervals
%   y0     initial condition (vector of length m)
%   f      function handle defining the problem
%            f = @(t,y) ....
%          where t = time, y = vector of length m
%          and output is a vector of length m
%
% Outputs:
%  y       iterate history (m by n+1 matrix)
function y = rk4(a,b,n,y0,f)
 y=zeros(length(y0),n+1);
 y(:,1) = y0;
 h=(b-a)/n;
 for i=1:n,
  ti=a+(i-1)*h;
  F1 = h*f(ti,y(:,i));
  F2 = h*f(ti+h/2,y(:,i)+F1/2);
  F3 = h*f(ti+h/2,y(:,i)+F2/2);
  F4 = h*f(ti+h  ,y(:,i)+F3);
  y(:,i+1) = y(:,i) + (F1 + 2*F2 + 2*F3 + F4)/6;
 end;
