function [T,Y] = rk4V(Fsys,Tspan,Y0,n)
%   [T,Y]=rkV4(F,[a b],Y0,n)
% Fixed step h=(b-a)/n Runge-Kutta for an ODEsystem 
% Y' = Fsys(t,Y), Y(a)) = Y0.
% Input:
% Fsys: Funtion that defines the system (as in EulerV)
% Tspan: Tinme interval [t0,t].
% Y0:    Vector of initial values at time point t0.
% n:     Nr. of steps
% Results:
% T:  Time column vector (T=linspace(a,b,n)';)
% Y:  Matrix, whose j^{th} column gives solution approximations, ...
%     at time point T_j
%     i^{th} row of [T Y]: [T(i)  Y(i,1), Y(i,2), ..., Y(i,N)]
%
% Example: Fsys=@(t,y)[y(1)-y(1)*y(2)/2;-2*y(2)+y(1)*y(2)]; 
%      y0=[20;10];
%      [T,Y]=rk4V(Fsys,[0,1],y0,20); plot(T,Y); Solutions' time-plots
%      figure; plot(Y(:,1),Y(:,2))  % Phase plane (2x2 autonomous system)
     
a=Tspan(1);b=Tspan(2);
h=(b-a)/n;
y0=[20;10];
N=length(Y0);    % Size of system 
Y=zeros(n+1,N);  % i^{th} row:  [Y(i,1), Y(i,2), ..., Y(i,N)]
T=a:h:b;
Y(1,:)=Y0';

for j=1:n,
  tj = T(j);  
  Yj = Y(j,:); 
  K1 = h*Fsys(tj,Yj)';
  K2 = h*Fsys(tj+h/2,Yj+K1/2)';
  K3 = h*Fsys(tj+h/2,Yj+K2/2)';
  K4 = h*Fsys(tj+h,Yj+K3)';
  Y(j+1,:) = Yj + (K1 + 2*K2 + 2*K3 + K4)/6;
end;
