function [tn,un]=eulero_esp_sys(odefun,tspan,y0,Nh)
% Eulero esplicito per sistemi 
% pbl di Cauchy del primo ordine
% [tn,un]=eulero_esp_sys(odefun,tspan,y0,Nh)
% Input: odefun = function handle di f(t,y) (restituisce un vettore
%                 colonna con numero di elementi pari al numero di eqz)
%        tspan = [t0, T] vettore con istante iniziale e finale
%        y0 = vettore delle condizioni iniziali del p. di Cauchy, 
%              (ha dimensione pari al numero di equazioni) 
%        Nh = numero di passi temporali
% Output: tn = vettore colonna dei tempi
%         un = matrice della soluzione numerica, ogni colonna restituisce
%             l'approssimazione della corrispondente componente 
%             della soluzione


% costruisco il vettore dei tempi (colonna)
tn=linspace(tspan(1),tspan(2),Nh+1)';
% condizione iniziale
y0=y0(:); % converte in vettore colonna un vettore y0
d=length(y0); % numero di equazioni del sistema
% inizializzo un
un=zeros(Nh+1,d);
% calcolo h=(T-t0)/Nh
h=(tspan(2)-tspan(1))/Nh;

% salvo la condizione iniziale in un
un(1,:)=y0.'; % trasposto valido anche per valori complessi

% ciclo in in
for n=1:Nh
    wn=un(n,:).';
    % applico un passo di Eulero esplicito
    w=wn+h*odefun(tn(n),wn);
    % salvo w all'interno di un
    un(n+1,:)=w.';
end
