function [t,y] = Euler(f,t0,T,y0,N)
%  [ts,ys] = Euler(f,t0,T,y0,N)
%
% Use N steps of Euler method to solve initial value problem
%  y'(t) = f(t,y) , y(t0) = y0
%
% INPUT:
%
% f(t,y)  is function on right hand side of differential equation
%           for inline function use Euler(f,...)
%           for m-file          use Euler('f',...)
%     t0  is initial time
%      T  is final time
%     y0  is initial value for y at time t0
%      N  is number of Euler steps
%
% OUTPUT:
%
%     ts  is vector of t-values
%     ys  is matrix of y-values: row i contains y-values at time t(i)
%
% EXAMPLE FOR ONE ODE:
%     y' = -y + t,  y(1)=2, find y(t) for t in [1,5] using 20 steps:
%
%     f = inline('-y + t', 't','y')
%     [ts,ys] = Euler(f, 1, 5, 2, 20);
%     [ts,ys]              % show table of t and y values
%     plot(ts,ys,'o-')     % plot solution
%
% EXAMPLE FOR SYSTEM OF ODEs:
%     y1' = -2*y1 + y2 + t
%     y2' = y1 - 3*y2, 
%     y1(1) = -1
%     y2(1) = 2 
% Find y(t) for t in [1,5] using 20 steps:
%
%     f = inline('[ -2*y(1) + y(2) + t; y(1) - 3*y(2) ]', 't','y')
%     [ts,ys] = Euler(f, 1, 5, [-1;2], 20);
%     [ts,ys]              % show table with columns for t, y1, y2
%     plot(ts,ys,'o-')     % plot solutions y1 and y2
%     legend('y1','y2')

t = zeros(N+1,1);              % initialize t and y with zero arrays
y = zeros(N+1,length(y0));

h = (T - t0)/N;
t(1) = t0;                     % initialize first rows of t,y with initial values
y(1,:) = y0.';

for i=1:N
  t(i+1) = t(i) + h;
  s = feval(f,t(i),y(i,:));    % evaluate s = f(t(i),y(i))
  y(i+1,:) = y(i,:) + h*s.';
end
