function [t,y] = eulsys(dydt,tspan,y0,h,varargin)
%
% eulsys: Euler ODE solver for systems of ODEs 
%   [t,y] = eulsys(dydt,tspan,y0,h,p1,p2,...):
%           uses Euler's method to integrate a system of ODEs
%
% input:
%   dydt = name of the M-file that evaluates the rhs of the ODEs 
%   tspan = [ti, tf] where ti and tf = initial and
%           final values of independent variable
%   y0 = initial value of dependent variable
%   h = step size
%   p1,p2,... = additional parameters used by dydt
%
% output:
%   t = vector of independent variable
%   y = matrix of solutions for dependent variables
%

if nargin<4,error('at least 4 input arguments required'),end

ti = tspan(1);tf = tspan(2);
if ~(tf>ti),error('upper limit must be greater than lower'),end

t = ti:h:tf; 
n = length(t);

% if necessary, add an additional value of t
% so that range goes from t = ti to tf
if t(n)<tf  
  t(n+1) = tf;
  n = n+1;
end

% determine size of y0
dims = size(y0);
% make sure y0 is a column vector
if dims(2)>dims(1)
  y0 = y0';
end
jmax = max(dims);

y = y0; % initialize y

% implement Euler method
for i = 1:n-1  % loop over time
  rhs = dydt(t(i),y(:,i),varargin{:});
  for j=1:jmax  % loop over equations
    y(j,i+1) = y(j,i) + rhs(j)*(t(i+1)-t(i));
  end
end

