

function [T,Y]=EuleroNewton(t0,tmax,n,y0,f,Jac)

%______________________________________________________________________________
%
% Funzione EuleroNewton
% ______________________
%
% Chiamata:   [T,Y]=EuleroNewton(t0,tmax,n,y0,f,Jac);
%
% Commento:
%
% La function EuleroNewton utilizza il metodo di Eulero del 
% primo ordine implicito per determinare i valori approssimati 
% della soluzione di un problema di Cauchy, nei nodi 
% t(i), i=1,2,...,n+1 contenuti in  un intervallo [t0,tmax].
% La function utilizza il metodo di Newton per la risoluzione
% del sistema non lineare ottenuto dalla discretizzazione.
% Il problema (scalare o vettoriale) e' del tipo:
% 
% 
% y'(t)=f(t,y(t))     con t in (t0,tmax]
% y(t0)=yo
%
%   
% dove:
% y(t)=(y1(t),y2(t),...,ym(t))
% f(t,y(t))=(f1(t,y(t)),f2(t,y(t)),...,fm(t,y(t))
% y(t0)=(y1(t0),y2(t0),...,ym(t0))=(yo1,yo2,...,yom)=yo
%
% con m=1 nel caso scalare e m>1 nel caso vettoriale.
%___________________________________________________________________________ 
%
%
% Parametri di input della function:
% ___________________________________
%
% t0,tmax  estremi dell'intervallo di integrazione.
%
%  n       numero di sottointervalli [t(i),t(i+1)], i=1,...,n
%          in cui si desidera suddividere l'intervallo [t0,tmax].    
%
% y0       vettore riga di lunghezza m contenente la condizione iniziale
%          y(t0)=yo=(yo1,yo2,...,yom).
%
% f        vettore colonna di dimensione m  in cui ogni componente e' 
%          una stringa contenente le espressioni f1, f2,..,fm scritte in
%          funzione delle variabili t e y; si puo' creare il vettore f 
%          utilizzando il comando 'strvcat', nel seguente modo: 
%          f=strvcat(f1,f2,...,fm). 
%
% Jac      Vettore colonna di dimensione m^2 in cui ogni componente 
%          e' una stringa contenente le espressioni delle derivate 
%          parziali delle fi, i=1,...,m, che costituiscono la matrice 
%          jacobiana; per m=2 Jac=[J11;J12;J21;J22] dove Jik, i,k=1,2,
%          e' la stringa della derivata di fi rispetto a yk;
%          Per creare il vettore Jac si puo' usare il comando 'strvcat', 
%          nel seguente modo: Jac=strvcat(J11,J12,J21,J22).
%          Si noti che per m=1 Jac e' una stringa contenente l'espressione
%          della derivata di f(t,y), fatta rispetto a y.

%
% Parametri di output della function:
% ___________________________________
%
% T        vettore colonna di dimensione n+1 contenente i nodi t0=t(1),
%          t(2),...,t(n),t(n+1)=tmax.
%
% Y        matrice di dimensione (n+1)x m la cui i-esima riga di lunghezza 
%          m contiene l'approssimazione della funzione y(t)=[y1(t),...,ym(t)], 
%          calcolata nel nodo t(i): y(t(i))=[y1(t(i)),y2(t(i)),...,ym(t(i))].
% 
%___________________________________________________________________________
%
%
% 1° esempio  di utilizzo della function  applicata al seguente 
% problema di Cauchy (caso scalare con m=1):
%
% y'(t)=-y+2cos(t)    t e' in (0,pi]
% y(0)=1.
%
% La soluzione vera del problema e':
% y(t)=sin(t)+cos(t).
%
%
%      clear all
%      t0=0; tmax=pi;
%      y0=1;
%      n=300;
%      f='-y+2*cos(t)';
%      Jac='-1'; %derivata della f fatta rispetto a y;
% 
%      [T,Y]=EuleroNewton(t0,tmax,n,y0,f,Jac);
%       
%___________________________________________________________________________
%
%
% 2° esempio  di utilizzo della function  applicata al seguente 
% problema di Cauchy (da ricondursi a un problema vettoriale del primo 
% ordine con m=2):
%
% 
% z''(t)=(32+2t^3-z(t)z'(t))/8    t e' in (1,3]   (*)
% z(1)=17
% z'(1)=-14.
%
% Il problema di Cauchy (*) e' del secondo ordine; per poter applicare
% la function EuleroNewton a tale problema e' necessario trasformare (*)
% in un problema del 1° ordine. 
% Ponendo z(t)=y1(t), il problema (*) si trasforma nel problema 
% equivalente:
%
% y1'(t)=y2(t)                    t e' in (1,3]   (**)
% y2'(t)=(32+2t^3-y1(t)y2(t))/8
% y1(1)=17
% y2(1)=-14.
%
% La soluzione vera di (**) e':
% y1(t)=t^2+16/t
% y2(t)=2t-16/t^2.
%
% La soluzione vera di (*) e' allora:
% z(t)=y1(t)=t^2+16/t.
% La y2(t) corrisponde alla derivata z'(t).
% 
% Il main relativo al problema dell'esempio potrebbe essere cosi' 
% configurato:
%
% t0=1; tmax=3;
% y0=[17 -14];
% n=200;
% f1='y(2)';
% f2='(32+2*t.^3-y(1)*y(2))./8';
% f=strvcat(f1,f2);
% Jac=strvcat('0','1','-y(2)/8','-y(1)/8'); %matrice jacobiana
% 
% [T,Y]=EuleroNewton(t0,tmax,n,y0,f,Jac);
% 
% % Osservare che in questo caso Y e' una matrice di dimensione 
% % n+1 x 2 contenente nella prima colonna i valori y1(t(i)), 
% % i=1,2,...,n+1, che sono  le approssimazioni nei nodi della 
% % soluzione del problema (*); la seconda colonna di Y contiene
% % i valori y2(t(i)) che sono le approssimazioni della derivata
% % prima della soluzione di (*).
%
%___________________________________________________________________________







h=(tmax-t0)/n;
m=length(y0);
Y=y0;
T=[t0:h:t0+n*h]';

toll=1.e-12;

for i=1:n
    Yold=Y(i,:);
    Y(i+1,:)=Yold+1;
    cont=0;y=Yold;t=T(i+1);
    
    while max(abs(Y(i+1,:)-Yold))>toll & cont<200
         for j=1:m
             k(j)=eval(f(j,:));
         end
         F=Yold-Y(i,:)-h*k;
         for k=1:m
             for j=1:m
                 JF(k,j)=-h*eval(Jac((k-1)*m+j,:));  
             end
         end
         JF=JF+eye(m);
         if abs(det(JF))<=eps
            disp('Matrice jacobiana con determinante nullo')
            return
         end
         diff=JF\(-F'); 
         Y(i+1,:)=Yold+diff';     
         y=Y(i+1,:);Yold=Y(i+1,:);
         cont=cont+1;
     end
     if cont>=200
        disp('Superato il numero massimo di iterazioni')
        return
     end
end

