function [t,y]=ForwardEuler(f,tspan,y0,n);
% FORWARDEULER solves system of ODEs using the Forward Euler method
%   [t,y]=ForwardEuler(f,tspan,y0,n) solves dy/dt=f(t,y) with initial
%   value y0 on the time interval tspan doing n steps of Forward Euler

dt=(tspan(2)-tspan(1))/n; 
t=tspan(1):dt:tspan(2);   
y(:,1)=y0(:);                           % colon to make column vector
for i=1:n,
  y(:,i+1)=y(:,i)+dt*f(t(i),y(:,i));
end;
t=t';y=y';                              % to return results in columns