

function [T,X,Y]=Eulero2(t0,tmax,n,x0,y0,f1,f2)

%_____________________________________________________________________________
%
% Funzione Eulero2
% ____________________
%
% Chiamata:   [T,X,Y]=Eulero2(t0,tmax,n,x0,y0,f1,f2);
%
% Commento:
%
% La function Eulero2 utilizza il metodo di Eulero in serie del
% primo ordine esplicito per determinare i valori approssimati
% della soluzione di un problema di Cauchy di due equazioni,
% nei nodi t(i), i=1,2,...,n+1 contenuti in  un intervallo 
% [t0,tmax]. Il problema e' del tipo:
%
%
% x'(t)=f1(t,x(t),y(t))     con t in (t0,tmax]  (*)
% y'(t)=f2(t,x(t),y(t))
% x(t0)=xo             
% y(t0)=yo.
%______________________________________________________________________
%
%
% 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].    
%
%  x0,y0    variabili reali contenenti le condizioni iniziali
%           x(t0)=xo e y(t0)=yo.
%
%  f1,f2    variabili stringa che contengono le espressioni (in 
%           funzione di t, x e y) di f1(t,x(t),y(t)) e di 
%           f2(t,x(t),y(t)) del problema (*).
%
%
% 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.
%
% X        vettore colonna di dimensione n+1, la cui i_esima componente 
%          contiene l'approssimazione  di x(t) calcolata nel nodo t(i).
%
% Y        vettore colonna di dimensione n+1, la cui i_esima componente 
%          contiene l'approssimazione di y(t) calcolata nel nodo t(i).
%
%________________________________________________________________________
%
%
% Esempio  di utilizzo della function  applicata al seguente 
% problema di Cauchy del primo ordine:  
% 
%
% x'(t)=x(t)*(r*y(t)-gamma)    t e' in (0,2]
% y'(t)=(-r*y(t)+gamma)*x(t)        
% x(0)=50
% y(0)=450.
%
%
% La soluzione vera del problema e':
% x(t)=(r*N-gamma)/(r+((r*N-gamma)/I0-r)*exp((-r*N+gamma)*t));
% y(t)=N-x(t)
%
% dove:
% r=0.02;gamma=0.03;
% N=500;I0=50;
%
% Osservare che si puo' utilizzare la function anche per risolvere
% numericamente il problema:
% z''(t)=f(t,z',z)
% z(t0)=zo
% z'(t0)=z1
%  
% dopo averlo opportunamente modificato in un problema vettoriale
% del primo ordine.
%
%
%
% clear all
% t0=0;tmax=2;
% x0=50;y0=450;
% n=200;
% f1='x.*(0.02*y-0.03)';
% f2='(-0.02*y+0.03).*x';
% 
% [T X Y]=Eulero2(t0,tmax,n,x0,y0,f1,f2);    
% 
%___________________________________________________________________________





h=(tmax-t0)/n;
T=t0;X=x0;Y=y0;
t=t0;x=x0;y=y0;

for i=1:n
  x=x+h*eval(f1);
  y=y+h*eval(f2);
  t=t+h;
  T=[T;t];
  X=[X;x];
  Y=[Y;y];
end
