%%
% Plots absolute-stability regions for ADAMS-BASHFORTH methods of orders 
% one through six. 
%
% Note: ADAMS-MOULTON demo code follows in a separate section below. 
% 
%

p = 6;

xmin = -2.5; xmax = .5; 
ymin = -1.5; ymax = 1.5;

n = 500;
h = (n+1)\2*pi;
theta = 0:h:2*pi;


figure(1); clf; axis square; hold on;
plot([xmin xmax],[0 0],'-k','LineWidth',1);
plot([0 0],[ymin ymax],'-k','LineWidth',1);
axis([xmin xmax ymin ymax]);grid on;axis equal;

ASRT = 'Absolute Stability Regions for Adams-Bashforth Methods';
title(ASRT,'FontSize',16);

for k = 1:p
    switch k
        case 1
% Forward Euler
plot( exp(i*theta)-1 ,'LineWidth',2);
title({ASRT;'Order 1 (Forward Euler)'});
pause;
        case 2
% 2nd order
plot( ((2\3)*exp(i*theta) - (2\1)).\(exp(i*(2*theta)) - exp(i*theta)),...
    'LineWidth',2);
title({ASRT;'Orders 1-2'});
pause;
        case 3 
% 3rd order
plot( ( (12\23)*exp(i*(2*theta)) - (3\4)*exp(i*theta) + 12\5).\...
    ( exp(i*(3*theta)) - exp(i*(2*theta)) ),'LineWidth',2);
title({ASRT;'Orders 1-3'});
pause;
        case 4
% 4th order
rat4 = @(theta)((24\55)*exp(i*(3*theta)) - (24\59)*exp(i*(2*theta)) + ...
    (24\37)*exp(i*theta) - 8\3).\(exp(i*(4*theta)) - exp(i*(3*theta)));

% Find indices of non-spurious points in the left-half plane. 
ind4 = find( (real(rat4(theta))<0)&(-.44<=imag(rat4(theta)))&...
    (imag(rat4(theta))<=.44) );

plot(rat4(theta(ind4)),'LineWidth',2);
title({ASRT;'Orders 1-4'});pause;

        case 5
% 5th order
rat5 = @(theta)( 720\( 1901*exp(i*(4*theta)) - 2774*exp(i*(3*theta)) + ...
    2616*exp(i*(2*theta)) -1272*exp(i*theta)+ 251 ) ).\ ...
    (exp(i*(5*theta)) - exp(i*(4*theta)));

% Find indices of non-spurious points in the left-half plane. 
ind5 = find( (real(rat5(theta))<0)&(-.23<=imag(rat5(theta)))&...
    (imag(rat5(theta))<=.23) );

plot(rat5(theta(ind5)),'LineWidth',2);
title({ASRT;'Orders 1-5'});pause;

        case 6
% 6th order
rat6 = @(theta)( 1440\( 4277*exp(i*(5*theta)) - 7923*exp(i*(4*theta)) + ...
    9982*exp(i*(3*theta)) - 7298*exp(i*(2*theta)) + 2877*exp(i*theta)- ...
    475 ) ).\ (exp(i*(6*theta)) - exp(i*(5*theta)));

% Find indices of non-spurious points in the left-half plane. 
ind6 = find( (real(rat6(theta))<0)&(-.12<=imag(rat6(theta)))&...
    (imag(rat6(theta))<=.12) );
plot(rat6(theta(ind6)),'LineWidth',2);
title({ASRT;'Orders 1-6'});

    end % End of switch. 
              
end % End of outer loop 

hold off;

%% 
% Plots absolute-stability regions for ADAMS-MOULTON methods of orders 
% one through six. 
% 
%

p = 6;

xmin = -7; xmax = 2.5; 
ymin = -4; ymax = 4;

n = 500;
h = (n+1)\2*pi;
theta = 0:h:2*pi;


figure(2); clf; axis square; hold on;
plot([xmin xmax],[0 0],'-k','LineWidth',1);
plot([0 0],[ymin ymax],'-k','LineWidth',1);
axis([xmin xmax ymin ymax]);grid on;axis equal;

ASRT = 'Absolute Stability Regions for Adams-Moulton Methods';
title(ASRT,'FontSize',16);

for k = 1:p
    switch k
        case 1
% Backward Euler
plot( exp(i*theta).\(exp(i*theta) - 1),'LineWidth',2);
title({ASRT;'Order 1 (Backward Euler, exterior of circle)'});
pause;
        case 2
% 2nd order
plot( ((2\1)*(exp(i*theta) +1)).\(exp(i*theta) - 1),'LineWidth',2);
title({ASRT;['Orders 1-2 (Backward Euler, exterior of circle;'...
    ' Trapezoidal Method, left half plane)']});
pause;
        case 3 
% 3rd order
plot( ( (12\1)*( 5*exp(i*(2*theta)) + 8*exp(i*theta) - 1) ).\...
    ( exp(i*(2*theta)) - exp(i*theta) ),'LineWidth',2);
title({ASRT;'Orders 1-3'});
pause;
        case 4
% 4th order
plot( ( (24\1)*( 9*exp(i*(3*theta)) + 19*exp(i*(2*theta)) - ...
    5*exp(i*theta) + 1) ).\( exp(i*(3*theta)) - exp(i*(2*theta)) ),...
    'LineWidth',2);
title({ASRT;'Orders 1-4'});
pause;
        case 5
% 5th order
plot( ( (720\1)*( 251*exp(i*(4*theta)) + 646*exp(i*(3*theta)) - ...
    264*exp(i*(2*theta)) + 106*exp(i*theta) - 19) ).\...
    ( exp(i*(4*theta)) - exp(i*(3*theta)) ),'LineWidth',2);
title({ASRT;'Orders 1-5'});
pause;
        case 6
%6th order
plot( ( (1440\1)*( 475*exp(i*(5*theta)) + 1427*exp(i*(4*theta)) - ...
    798*exp(i*(3*theta))+482*exp(i*(2*theta))-173*exp(i*theta)+27) ).\...
    ( exp(i*(5*theta)) - exp(i*(4*theta)) ),'LineWidth',2);
title({ASRT;'Orders 1-6'});
%pause;

    end % End of switch. 
              
end % End of outer loop 

hold off;
