

function [T,Y]=Eulsis(t0,tmax,n,y0,A,b)
%_____________________________________________________________________
%
% Funzione Eulsis
% ____________________
%
% Chiamata:   [T,Y]=Eulsis(t0,tmax,n,y0,A,b);
%
% Commento:
%
% La function Eulsis utilizza il metodo di Eulero del primo 
% ordine implicito per determinare i valori approssimati
% della soluzione di un problema di Cauchy a coefficienti 
% costanti, nei nodi t(i), i=1,2,...,n+1 contenuti in  un 
% intervallo [t0,tmax].
% Il problema e' del tipo:
%
%    
% y'(t)=Ay(t)+b(t)     con t in (t0,tmax]    (*)
% y(t0)=yo
%
%   
% dove:
% y(t)=(y1(t),y2(t),...,ym(t))
% A e' la matrice mxm dei coefficienti (costanti) del problema 
% b(t)=(b1(t), b2(t),..., bm(t))
% y(t0)=(y1(t0),y2(t0),...,ym(t0))=(yo1,yo2,...,yom)=yo
% con m>=1.
%______________________________________________________________________
%
%
% 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).
%
%  A        matrice di dimensione mxm dei coefficienti del problema.
%           
%  b        vettore colonna di dimensione m  in cui ogni componente e' 
%           una stringa contenente le espressioni b1, b2,..,bm dei termini
%           noti del problema, scritte in funzione della variabile t; 
%           si puo' creare il vettore b utilizzando il comando 'strvcat', 
%           nel seguente modo: b=strvcat(b1,b2,...,bm). 
%
%
%
% 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 del primo ordine:  
% 
%
% y1'(t)=y1(t)+9y2(t)    t e' in (0,1]
% y2'(t)=4y1(t)+y2(t)        
% y1(0)=3
% y2(0)=1
% 
% con soluzione vera:
% y1(t)=3/4(exp(-5t)+3exp(7t))
% y2(t)=1/2(-exp(-5t)+3exp(7t)).
%
%
%
% clear all
% t0=0;tmax=1;
% y0=[3 1];
% n=100;
% A=[1 9 ; 4 1];
% b=strvcat('0','0');
% 
% [T Y]=Eulsis(t0,tmax,n,y0,A,b);    
%
%________________________________________________________________________
%
% 2° esempio  di utilizzo della function  applicata al seguente 
% problema di Cauchy del primo ordine:  
% 
%
% y1'(t)=-2y1(t)+3y2(t)-t(t-5)                    t e' in (1,2]
% y2'(t)=  y1(t)-4y2(t)-exp(-2t)+3t^2-2t-1        
% y1(1)=1+exp(-2)
% y2(1)=0
% 
% con soluzione vera:
% y1(t)=exp(-2t)+t^2
% y2(t)=t(t-1).
%
%
%
% clear all
% t0=1;tmax=2;
% y0=[1+exp(-2) 0];
% n=100;
% A=[-2 3 ; 1 -4];
% b1='-t.*(t-5)';
% b2='-exp(-2*t)+3*t.^2-2*t-1';
% b=strvcat(b1,b2);
% 
% [T Y]=Eulsis(t0,tmax,n,y0,A,b);    
%
%___________________________________________________________________________








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

Y=y0;Yi=y0';
for i=1:n
    t=T(i+1);
    for j=1:m
        bi(j)=eval(b(j,:));
    end
    bsi=Yi+h*bi';
    Yi=Mat\bsi; 
    Y=[Y;Yi'];
end
