%
% 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   = matrix whose columsn 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] = Heun( FUN, t0, y0, h, N )
  % pre allcation of the output vector T and Y
  T = zeros(N+1,1) ;
  Y = zeros(size(y0,1),N+1) ; % size(y0,1) number of rows of y0

  % set T(1) and Y(1) with the initial condition
  T(1)   = t0 ;
  Y(:,1) = y0 ; %  Y(:,1) access the first column of matrix Y
  
  % loop performing the Euler Steps
  for k=1:N
    T(k+1)   = T(k) + h ;
    K1       = feval( FUN, T(k), Y(:,k) ) ;
    ETA      = Y(:,k) + h * K1 ; % Euler step y(k+1) = y(k) + h * f(t(k),y(k))
    K2       = feval( FUN, T(k+1), ETA ) ;
    Y(:,k+1) = Y(:,k) + (h/2) * (K1+K2) ; % Heun correction step 
                                          % y(k+1) = y(k) + (h/2) * (f(t(k),y(k))+f(t(k+1),y(k+1)))
  end

end
