%
% This routine approximate the solution of the IVP (Initial Value Problem)
%
%  /  y' = f(t,y)
%  |
%  |  y(a) = y[0]
%   \
%
%  PARAMETERS:
%   FUN = reference to user defined function
%   t0  = Starting integration time
%   y0  = Initial data of the IVP
%   h   = time step
%   N   = number of step to be performed
%
%  OUTPUT
%   T   = vector containing the nodes where 
%         the approximate solution is computed
%         t = t0, t0+h, t0+3*h,..... t0+N*h
%   Y   = vector containing the approximate solution at nodes T, i.e
%         Y(1) ~ y(t0), 
%         Y(2) ~ y(t0+h),
%         Y(3) ~ y(t0+2*h),
%         .....
%         Y(N+1) ~ y(t0+N*h)
%
function [T,Y] = Euler( FUN, t0, y0, h, N )
  % pre allcation of the output vector T and Y
  T = zeros(N+1,1) ;
  Y = zeros(N+1,1) ;
  
  % set T(1) and Y(1) with the initial condition
  T(1) = t0 ;
  Y(1) = y0 ;
  
  % loop performing the Euler Steps
  for k=1:N
    T(k+1) = T(k) + h ;
    Y(k+1) = Y(k) + h * feval( FUN, T(k), Y(k) ) ;
    % feval( FUN, T(k), Y(k) ) == FUN( T(k), Y(k) )
  end

end
