function [T,Y] = rk4V(fnimi,Tspan,Y0,n)
     %   [T,Y]=rkV4(F,[a b],Y0,n)
% Kiinteän askelpituuden h=(b-a)/n Runge-Kutta-ratkaisu systeemille:
% Syötteet:
% Y' = F(t,Y), Y(a)) = Y0.
% fnimi: Funktio, joka määrittelee systeemin (vrt. EulerV)
% Tspan: aikaväli [t0,t].
% Y0:    Alkuarvo(vektori) pisteessä t0.
% n:     ASkelten lukumäärä
% Tulokset:
% T:  Aika(pysty)vektori (T=linspace(a,b,n)';)
% Y:  Matriisi, jonka j:s sarake ilmaisee ratkaisufunktion
%           y_j arvot T-aikapisteissä.
%           i:s rivi: T(i)  Y(i,1), Y(i,2), ..., Y(i,N)
% Siis jos on 2x2-systeemi, niin Y:ssä on 2 saraketta.
%
% Esim: Fsys=inline('[y(1)-y(1)*y(2)/2;-2*y(2)+y(1)*y(2)]','t','y'); 
%       y0=[20;10];
%       [T,Y]=rk4V(Fsys,[0,1],y0,20); plot(T,Y); % aikakuva
%       figure; plot(Y(:,1),Y(:,2))  % faasikuva.
     %  Huom! Jos Fsys määritellään m-funktioksi, on kutsussa oltava 
     %  'Fsys' tai @Fsys.

a=Tspan(1);b=Tspan(2);
h=(b-a)/n;
N=length(Y0);    % Systeemin koko. 
Y=zeros(n+1,N);  % i:s rivi:  Y(i,1), Y(i,2), ..., Y(i,N)
T=a:h:b;
Y(1,:)=Y0';

for j=1:n,
  tj = T(j);   % lyhyemmät...
  Yj = Y(j,:); % .. nimet
  K1 = h*feval(fnimi,tj,Yj)';
  K2 = h*feval(fnimi,tj+h/2,Yj+K1/2)';
  K3 = h*feval(fnimi,tj+h/2,Yj+K2/2)';
  K4 = h*feval(fnimi,tj+h,Yj+K3)';
  Y(j+1,:) = Yj + (K1 + 2*K2 + 2*K3 + K4)/6;
end;
