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

%    modified from eulerode.m on June 9, 2017 by Whenson Zhou
%  using a hint from rkr4.m

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
% pre allocate the solution matrix to be an n rows, column depends on the number of column of y0
y = ones(n,1)*y0; 

% implement Euler's method
% add: for all y solution matrix, and transpose the dydt into a row matrix

for i = 1:n-1 
y(i+1,:) = y(i,:) + dydt(t(i),y(i,:),varargin{:})'*(t(i+1)-t(i));
end
