%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Forward Euler method for the solution of ODEs of the form:
% y'(t)  =  f(t,y(t))
% The entries of this function are:
% a) Initial time
% b) Final time
% c) Number of time steps between initial and final time steps
% d) The function f(t,y(t))
% e) Initial condition y(t=0)=y0
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
function [y,t]  =  ForwardEuler(t0,tf,NTime,f,DfDy,y0)

Dt  =  (tf-t0)/(NTime-1);  %  Time step constante
t   =  (linspace(t0,tf,NTime))';  %  Variable independiente, normalmente tiempo

NVariables  =  length(y0);
y           =  zeros(NTime,NVariables);             %  Inicialización de la variable dependiente

y(1,:)      =  y0;

for itime=2:length(t)
    y(itime,:)  =  y(itime-1,:) + Dt*f(t(itime-1),y(itime-1,:))';
    if nargin==6
       if size(y0,1)==1
          %----------------------------------------------------------------
          % Problema escalar
          %----------------------------------------------------------------
          JacobianF  =  DfDy(t(itime-1),y(itime-1,:));
          LimitingDt  =  2/abs(JacobianF);
          if Dt>LimitingDt
             warning('Dt larger than stability limit. Results may be unstable\n') 
          end
       else
          %----------------------------------------------------------------
          % Problema vectorial
          %----------------------------------------------------------------
          JacobianF   =  max(eig(DfDy(t(itime-1),y(itime-1,:))));
          LimitingDt  =  2/abs(JacobianF);
          if Dt>LimitingDt
             warning('Dt larger than stability limit. Results may be unstable\n') 
          end
       end

    end
end

 
