

function [T,Y]=Adamsbm(t0,tmax,n,y_inn,f)

%___________________________________________________________________________
%
%  Funzione Adamsbm
%  ________________
%
%
%  Chiamata: [T,Y]=Adamsbm(t0,tmax,n,y_inn,f);   
%
%
% La function Adamsbm utilizza il metodo di Adams-Bashforth-Moulton  
% del quarto ordine 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].
% 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].    
%
% y_inn    matrice di dimensione 4 x m la cui i-esima riga di 
%          lunghezza m contiene l'approssimazione ottenuta, ad 
%          esempio, dal metodo di Runge-Kutta del quarto ordine, 
%          della soluzione calcolata nel nodo t(i): 
%          y(ti)=[y1(ti),y2(ti),...,ym(ti)], i=0,1,2,3.
% 
% 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). 
%
%
% 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'=y(y-2)t                   t e' in (0,1]
% y(0)=1
%
% con soluzione: y=2/(exp(t^2)+1).
%
%
%
%
%      clear all
%      t0=0; tmax=1;
%      y0=1;
%      n=100;
%      f='y.*(y-2).*t';
%      h=(tmax-t0)/n;
%
%      [t_inn,y_inn]= Rungekutta4(t0,t0+3*h,3,y0,f);
%      [T,Y]=Adamsbm(t0,tmax,n,y_inn,f);
%
%
%___________________________________________________________________________
%
%
% 2° esempio  di utilizzo della function  applicata al seguente 
% problema di Cauchy (problema vettoriale del primo ordine con m=2):
%
%
% y1'(t)=y1(t)*(r*y2(t)-gamma)    t e' in (0,2]
% y2'(t)=(-r*y2(t)+gamma)*y1(t)        
% y1(0)=50
% y2(0)=450.
%
%
% La soluzione vera del problema e':
% y1(t)=(r*N-gamma)/(r+((r*N-gamma)/I0-r)*exp((-r*N+gamma)*t));
% y2(t)=N-y1(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;
% n=200;
% y0=[50 450];
% h=(tmax-t0)/n;
% f1='y(1).*(0.02*y(2)-0.03)';
% f2='(-0.02*y(2)+0.03).*y(1)';
% f=strvcat(f1,f2);
%
% [t_inn,y_inn]= Rungekutta4(t0,t0+3*h,3,y0,f);
% [T,Y]=Adamsbm(t0,tmax,n,y_inn,f);
%
% % 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 
% % prima componente della soluzione del problema (*); la seconda
% % colonna di Y contiene i valori y2(t(i)) che sono le 
% % approssimazioni della seconda componente della soluzione di (*).
%
%___________________________________________________________________________






        

h=(tmax-t0)/n;
[l,m]=size(y_inn);      
v=zeros(1,m);
T=[t0:h:t0+n*h]';
Y=y_inn;
r=4;

for i=1:r
    t=T(i);y=Y(i,:);
	for j=1:m
   	    fun(i,j)=eval(f(j,:));
    end
end
fun=[fun;v];

coef1=[-9 37 -59 55 0]/24;
coef2=[0 1 -5 19 9]/24;
for i=r:n
    p=Y(i,:)+h*coef1*fun;                    %predictor
    
    t=T(i+1);y=p;                            
    for j=1:m
   	    fun(r+1,j)=eval(f(j,:));             %evaluation
    end
    Y(i+1,:)=Y(i,:)+h*coef2*fun;             %corrector
    
    y=Y(i+1,:);                              
    for j=1:m
   	    fun(r+1,j)=eval(f(j,:));             %evaluation
    end
    fun=[fun(2,:);fun(3,:);fun(4,:);fun(5,:);v];
end    

