function [t,x]=euler_forw(odefun,tspan,x0,h,varargin)
% [t,x]=FEULER(odefun,tspan,x0,Nh) 
%   Risolve l'IVP
%   x'(t)=odefun(t, x(t)), t in [tspan(1), tspan(2)]
%   tspan(1)=x0
%
%   usando il metodo di Eulero in avanti su una griglia 
%   equispaziata con intervalli di ampiezza h.
%   La funzione ODEFUN(T,X) ritorna una matrice, la cui generica riga
%   e' l'approssimazione della soluzione in un punto t_i
%   contenente f(t,x), della stessa dimensione di x.

%   [t,x] = FEULER(ODEFUN,TSPAN,X0,H,P1,P2,...) passa
%   i parametri addizionali P1,P2,... alla funzione
%   ODEFUN, definita come ODEFUN(T,X,P1,P2...).

n = fix((tspan(2)-tspan(1))/h) + 1;
t = linspace(tspan(1),tspan(1)+(n-1)*h,n)';

w=x0(:); x=w.'; % w e' un vettore colonna

for i = 1:n-1
 w=w+h*feval(odefun,t(i),w,varargin{:});
 x = [x; w.'];
end

return
