function [tn,un]=eulero_imp_sys(odefun,tspan,y0,Nh,varargin)
% Eulero implicito per sistemi 
% pbl di Cauchy del primo ordine
% [tn,un]=eulero_imp_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


% controllo il numero di input
if nargin==4 % non ho passato parametri opzionali
    % li definisco io
    tol=1e-8; kmax=20;
elseif nargin==5
    tol=varargin{1};
    kmax=20;
elseif nargin==6
    tol=varargin{1};
    kmax=varargin{2};
end

% 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

% ad ogni iterazione di Eulero devo risolvere un 
% sistema di eqz non lineari, utilizzo Broyden
B0=eye(d); % matrice iniziale per Broyden
for n=1:Nh
    wn=un(n,:).';
    % definisco la funzione non lineare
    r=@(x)x-wn-h*odefun(tn(n+1),x);
    % calcolo lo zero con Broyden
    [w,res,it]=broyden(r,B0,wn,tol,kmax);
    un(n+1,:)=w.';
end





