%Computation of the A-Stability Region of explicit Adams-Moulton metohds
clear, clf reset

% 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 implicito
beta(1,1)=1; 
% metodo dei trapezi
beta(2,1:2)=[1, 1]/2;
% metodo a 2 passi
beta(3,1:3)=[5, 8,-1]/12;
% metodo a 3 passi
beta(4,1:4)=[9, 19, -5, 1]/24;
% metodo a 4 passi
beta(5,1:5)=[251, 646, -264, 106, -19]/720;


% La frontiera della regione di stabilita' assoluta viene disegnata
% disegnando i punti z=rho(om)/sigma(om), al variare di om sulla
% circonferenza unitaria, dove rho(.) e sigma(.) sono i 2 polinomi
% associati al LMM
r=[0:200];
om=cos(2*pi/200*r)+i*sin(2*pi/200*r);
z=(om-1)./(polyval([beta(1,1) 0],om));
plot(z,'r.')
hold on
z=(om-1)./(polyval([beta(2,1:2)],om));
iz=(abs(z)<30);
plot(z(iz),'g.')
z=(om.^2-om)./(polyval([beta(3,1:3)],om));
plot(z,'b.')
z=(om.^3-om.^2)./(polyval([beta(4,1:4)],om));
plot(z,'k.')
z=(om.^4-om.^3)./(polyval([beta(5,1:5)],om));
plot(z,'y.')


grid on
axis equal
plot([-2.5 2],[0 0 ],'-k')
plot([0 0],[-2 2],'-k')

legend('Eulero implicito','Trapezi','k=2','k=3','k=4')
title('Metodi impliciti di Adams-Moulton, Stability Regions')
