function [T,Y]=eulerV(Fsys,Tspan,ya,n)
% call:   [T,Y]=eulerV(Fsys,Tspan,ya,n);
%  or:    [T,Y]=eulerV(@Fsys,Tspan,ya,n); % if Fsys is an m-file-function.
% Inputs:
% - Fsys : Function handle or f(t,y) defining the system y'=f(t,y).
%          Same as for ode45  (help ode45)
% - Tspan:   is of the form [a b]  
% - ya:      Vector of initial value (y(a)=ya) (row or column)
% - n:       Nr. of subintervals of [a,b] 
%
% Ouputs: 
%  - T:   Time points, column vector
%  - Y:   Matrix whose jth column has the values of
%         solution approximations y_j at T-timepoints.
%   (i:th row of [T Y]: [T(i) Y(i,1), Y(i,2), ..., Y(i,N)] )
%
% Example: 
%       N=20  % 200, ...
%       tf=6;
%       rabfox=@(t,y,alpha)[2*y(1)-alpha*y(1)*y(2);-y(2)+alpha*y(1)*y(2)]; 
%       y0=[300;150];  
%       alpha=0.01;
%       [T,Y]=eulerV(@(t,y)rabfox(t,y,alpha),[0,tf],y0,N);   
%       plot(T,Y)  % Solutions in (t,y)-axes
%       figure
%       plot(Y(:,1),Y(:,2)) % Phase plane trajectory (2x2 autonomous sys.)
ya=ya(:)';   % Make row vector 
a=Tspan(1);b=Tspan(2);
h=(b-a)/n;
N=length(ya);
Y=zeros(n+1,N);  % j^{th} col:  Y(1,j), Y(2,j), ..., Y(N,j)
T=a:h:b;
Y(1,:)=ya;       % First row

for i=1:n
   Y(i+1,:)=Y(i,:)+h*(Fsys(T(i),Y(i,:))');                                              
end;
