function [t, x]= adams(fname, tspan, x0, h, k)
% Approssima una soluzione di x'=f(t,x), x(t0)=x0
% nell'intervallo [t0,T].
% con il metodo di Adams-Bashforth esplicito a a k passi,
% per k=1,...,6: x(i+1)=x(i) + h sum_{j=1:k} beta_j f_{i-j+1}
% Il metodo a k passi e' convergente di ordine k
%
% Input:
% fname: nome della funzione f(t,x):RxR^n ----> R^n
% t0, x0: valore iniziale, x0=x(t0)
% T: estremo destro intervallo d'integrazione
% h: lunghezza del passo
% k: numero di passi del metodo
%
% Output:
% t, x: array tali che x(i) e' un'approssimazione di x(t(i)) e
% t(i) = t0 + (i-1)h, i=1:n+1
%
% Le prime k approssimazioni sono calcolate con il metodo di Adams
% di ordine 1, 2,..., k-1
%
% Si provi la funzione f5.m
t0=tspan(1);
T=tspan(2);
t=[t0:h:T];
n=length(t);
x=zeros(1,n);
f=zeros(1,n); 

% la matrice beta contiene, sulla riga k, i coefficienti beta_i del metodo
% a k passi, completati con zeri fino a 6
beta=zeros(6);
% metodo di eulero
beta(1,1)=1; 
% metodo a 2 passi
beta(2,1:2)=[3, -1]/2;
% metodo a 3 passi
beta(3,1:3)=[23, -16, 5]/12;
% metodo a 4 passi
beta(4,1:4)=[55, -59, 37, -9]/24;
% metodo a 5 passi
beta(5,1:5)=[1901, -2774, 2616, -1274, 251]/720;
% metodo a 6 passi
beta(6,1:6)=[4277, -7923, 9982, -7298, 2877, -475]/1440;

% Calcola le prime k approssimazioni
x(1)=x0;
f(1)=feval(fname, t(1), x(1));
for m=2:k
    x(m)=x(m-1)+h* beta(m-1, 1:m-1)*f(m-1:-1:1)';
    f(m)=feval(fname, t(m), x(m));
end
% calcola le successive approssimazioni
for m=k+1:n
    x(m)=x(m-1)+h* beta(k, 1:k)*f(m-1:-1:m-k)';
    f(m)=feval(fname, t(m), x(m));
end

