function [T,Z] = rks4(Fn,a,b,Za,m)
% [T,Y] = rks4(f,a,b,ya,m)
% Runge-Kutta solution for the system
% of DE's   Z' = F(t,Z) with Z(a) = Za.
% f  is the function, input.
% a  is the left endpoint, input.
% b  is the right endpoint, input.
% Za is the initial condition vector, input.
% m  is the number of steps, input.
% T  is the parameter vector, output.
% Z  is the matrix of the vector solutions, output.
h = (b - a)/m;
T = zeros(1,m+1);
Z = zeros(m+1,length(Za));
T(1) = a;
Z(1,:) = Za;
for j=1:m,
  tj = T(j);
  Zj = Z(j,:);
  K1 = h*feval(Fn,tj,Zj);
  K2 = h*feval(Fn,tj+h/2,Zj+K1/2);
  K3 = h*feval(Fn,tj+h/2,Zj+K2/2);
  K4 = h*feval(Fn,tj+h,Zj+K3);
  Z(j+1,:) = Zj + (K1 + 2*K2 + 2*K3 + K4)/6;
  T(j+1) = a + h*j;
end
