function [t,y] = eulersys(t0,y0,t_end,h,fcn)
%
% function [t,y]=eulersys(t0,y0,t_end,h,fcn)
%
% Solve the initial value problem of a system
% of first order equations
%    y' = f(t,y),  t0 <= t <= b,  y(t0)=y0
% Use Euler's method with a stepsize of h.  
% The user must supply a program to compute the
% right hand side function with some name, say
% deriv, and a first line of the form
%   function ans=deriv(t,y)
% A sample call would be
%   [t,z]=eulersys(t0,z0,b,delta,'deriv')
%
% The program automatically determines the
% number of equations from the dimension of
% the initial value vector y0.
%
% Output:
% The routine eulersys will return a vector t
% and matrix y. The vector t will contain the
% node points in [t0,t_end]:
%   t(1)=t0, t(j)=t0+(j-1)*h, j=1,2,...,N
% The matrix y is of size N by m, with m the
% number of equations. The i-th row y(i,:) will
% contain the estimates of the solution Y 
% at the node points in t(i).
%
m = length(y0);
n = fix((t_end-t0)/h)+1;
t = linspace(t0,t0+(n-1)*h,n)';
y = zeros(n,m);
y(1,:) = y0;
for i = 2:n
  y(i,:) = y(i-1,:) + h*feval(fcn,t(i-1),y(i-1,:));
end

end % eulersys