%Computation of the A-Stability Region of implicit  BDF metohds
clear, clf reset

% la matrice alfa contiene, sulla riga k, i coefficienti alfa_i del metodo
% a k passi, completati con zeri fino a 6
alfa=zeros(6);
% metodo di eulero implicito
alfa(1,1:2)=[1 -1]; 
% metodo a 2 passi
alfa(2,1:3)=[3, -4, 1]/3;
% metodo a 2 passi
alfa(3,1:4)=[11, -18, 9, -2]/11;
% metodo a 3 passi
alfa(4,1:5)=[25, -48, 36, -16, 3]/25;
% metodo a 4 passi
alfa(5,1:6)=[137, -300, 300, -200, 75, 12]/137;


% 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=polyval([alfa(1,1:2)],om)./(om);
plot(z,'r.')
hold on
z=polyval([alfa(2,1:3)],om)./(2/3*om.^2);
plot(z,'g.')
z=polyval([alfa(3,1:4)],om)./(6/11*om.^3);
plot(z,'b.')
z=polyval([alfa(4,1:5)],om)./(12/25*om.^4);
plot(z,'k.')
z=polyval([alfa(5,1:6)],om)./(60/147*om.^5);
plot(z,'y.')


grid on
axis equal
plot([-5 20],[0 0 ],'-k')
plot([0 0],[-10 10],'-k')

legend('Eulero implicito','k=2','k=3','k=4','k=5')
title('Metodi impliciti BDF, Stability Regions')
