function [tn,un]=ee_adattivo_sys(odefun,tspan,y0,tol,hmin)
% Eulero esplicito adattivo per sistemi
% Input: odefun=@(t,y): function handle della funzione del sistema, deve
%        restiture un vettore colonna di d componenti, tante quante sono le 
%        equazioni del sistema
%        tspan=[t0,T]: vettore con gli estremi dell'intervallo di
%        simulazione
%        y0= vettore colonna delle condizioni iniziali. E' un vettore di di
%        componente, d e' la dimensione del sistema
%        tol= tolleranza per lo stimatore dell'errore
%        hmin= minimo valore di h accettato
% Output: tn= vettore colonna dei tempi
%         un= matrice con di colonne, ognuna contenente una componente
%         della soluzione del sistema
%   [tn,un]=ee_adattivo_sys(odefun,tspan,y0,tol,hmin)

d=length(y0);
y0=y0(:);
h=diff(tspan)/10; % h iniziale
tn=[tspan(1)];
un=[y0.']; 
% definisco le variabili t e u 
% che contengono l'ultimo tempo e l'ultima soluzione
t=tn(1); u=un(1,:).';
while t< tspan(2)
    % soluzione di EE con 1 passo h
    w=u+h*odefun(t,u);
    % soluzione di EE con 2 passi h/2
    wbar=u+h/2*odefun(t,u);
    wtilde=wbar+h/2*odefun(t+h/2,wbar);
    % controllo sullo stimatore dell'errore
    if norm(w-wtilde,"inf")<= 0.5*norm(un,inf)*tol || h< hmin
        % accetto h e la soluzione
        un=[un;wtilde.'];
        tn=[tn;t+h];
        h=h*2; % h per il prossimo passo
        % aggiorno t e u per il passo successivo
        t=tn(end); u=un(end,:).';
    else
        h=h/2;
    end
end