% Plot stability regions of various linear multistep methods.

 T = exp(linspace(0,2i*pi,500));
 fillcolor = .71*[1 1 1];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Forward Euler: stability is interior of boundary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 alpha = [1 -1];
 beta  = [0 1];
 ax = 2.5*[-1 1 -1 1];
 figure(1), clf
 stabcurve = polyval(alpha,T)./polyval(beta,T);
 obj = fill(real(stabcurve),imag(stabcurve),fillcolor); hold on
 set(obj,'edgecolor',fillcolor)
 plot(ax([1 2]),[0 0],'k-','linewidth',.5); hold on
 plot([0 0],ax([3 4]),'k-','linewidth',.5); 
 axis equal, axis(ax)
 set(gca,'fontsize',20)
 title('Forward Euler Method', 'fontsize',20)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Backward Euler: stability is exterior of boundary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 alpha = [1 -1];
 beta  = [1 0];
 ax = 2.5*[-1 1 -1 1];
 figure(2), clf
 axis equal, axis(ax)
 fill(ax([1 2 2 1]),ax([3 3 4 4]),fillcolor); hold on
 stabcurve = polyval(alpha,T)./polyval(beta,T);
 obj = fill(real(stabcurve),imag(stabcurve),[1 1 1]);
 axis equal, axis(ax);
 set(obj,'edgecolor',fillcolor)
 plot(ax([1 2]),[0 0],'k-','linewidth',.5); hold on
 plot([0 0],ax([3 4]),'k-','linewidth',.5); 
 set(gca,'fontsize',20)
 title('Backward Euler Method', 'fontsize',20)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2nd Order Adams Bashforth: stability is interior of boundary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 alpha = [1 -1 0];
 beta  = [0 1.5 -.5];
 ax = 2.5*[-1 1 -1 1];
 figure(3), clf
 stabcurve = polyval(alpha,T)./polyval(beta,T);
 obj = fill(real(stabcurve),imag(stabcurve),fillcolor); hold on
 set(obj,'edgecolor',fillcolor)
 plot(ax([1 2]),[0 0],'k-','linewidth',.5); hold on
 plot([0 0],ax([3 4]),'k-','linewidth',.5); 
 axis equal, axis(ax)
 set(gca,'fontsize',20)
 title('2nd Order Adams-Bashforth Method', 'fontsize',20)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2nd Order Adams Moulton: stability is exterior of boundary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 alpha = [1 -1];
 beta  = [.5 .5];
 ax = 2.5*[-1 1 -1 1];
 figure(4), clf
 axis equal, axis(ax);
 fill(ax([1 2 2 1]),ax([3 3 4 4]),fillcolor); hold on
% stabcurve = polyval(alpha,T)./polyval(beta,T);
 stabcurve = [ax(3)*1i ax(4)*1i ax(4)*1i+ax(2) ax(3)*1i+ax(2)];
 obj = fill(real(stabcurve),imag(stabcurve),[1 1 1]);
 set(obj,'edgecolor',fillcolor)
 plot(ax([1 2]),[0 0],'k-','linewidth',.5); hold on
 plot([0 0],ax([3 4]),'k-','linewidth',.5); 
 axis equal, axis(ax);
 set(gca,'fontsize',20)
 title('2nd Order Adams-Moulton (Trapezoid) Method', 'fontsize',20)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4th Order Adams Bashforth: stability is interior of boundary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 alpha = [1 -1 0 0 0];
 beta  = [55 -59 37 -9]/24;
 ax    = [-4 1 -2.5 2.5];
 figure(5), clf
 stabcurve = polyval(alpha,T)./polyval(beta,T);
 stabcurve = stabcurve(real(stabcurve)<=0);
 obj = fill(real(stabcurve),imag(stabcurve),fillcolor); hold on
 set(obj,'edgecolor',fillcolor)
 plot(ax([1 2]),[0 0],'k-','linewidth',.5); hold on
 plot([0 0],ax([3 4]),'k-','linewidth',.5); 
 axis equal, axis(ax)
 set(gca,'fontsize',20)
 title('4th Order Adams-Bashforth Method', 'fontsize',20)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4th Order Adams Moulton: stability is interior of boundary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 alpha = [1 -1 0 0];
 beta  = [9 19 -5 1]/24;
 ax = [-4 1 -2.5 2.5];
 figure(6), clf
 stabcurve = polyval(alpha,T)./polyval(beta,T);
 obj = fill(real(stabcurve),imag(stabcurve),fillcolor); hold on
 set(obj,'edgecolor',fillcolor)
 plot(ax([1 2]),[0 0],'k-','linewidth',.5); hold on
 plot([0 0],ax([3 4]),'k-','linewidth',.5); 
 axis equal, axis(ax)
 set(gca,'fontsize',20)
 title('4th Order Adams-Moulton Method', 'fontsize',20)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1-step BDF: stability is exterior of boundary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 alpha = [1 -1];
 beta  = [.5 .5];
 ax = [-4 12 -8 8];
 figure(7), clf
 plot(ax([1 2]),[0 0],'k-','linewidth',.5); hold on
 plot([0 0],ax([3 4]),'k-','linewidth',.5);
 axis equal, axis(ax);
 fill(ax([1 2 2 1]),ax([3 3 4 4]),fillcolor); hold on
% stabcurve = polyval(alpha,T)./polyval(beta,T);
 stabcurve = [ax(3)*1i ax(4)*1i ax(4)*1i+ax(2) ax(3)*1i+ax(2)];
 obj = fill(real(stabcurve),imag(stabcurve),[1 1 1]);
 axis equal, axis(ax);
 set(obj,'edgecolor',fillcolor)
 plot(ax([1 2]),[0 0],'k-','linewidth',.5); hold on
 plot([0 0],ax([3 4]),'k-','linewidth',.5); 
 set(gca,'fontsize',20)
 title('1-step Backward Difference Formula (Trapezoid)', 'fontsize',20)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2-step Backward Difference Formula: stability is exterior of boundary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 alpha = [3 -4 1];
 beta  = [2 0 0];
 ax = [-4 12 -8 8];
 figure(8), clf
 axis equal, axis(ax)
 fill(ax([1 2 2 1]),ax([3 3 4 4]),fillcolor); hold on
 stabcurve = polyval(alpha,T)./polyval(beta,T);
 obj = fill(real(stabcurve),imag(stabcurve),[1 1 1]);
 axis equal, axis(ax);
 set(obj,'edgecolor',fillcolor)
 plot(ax([1 2]),[0 0],'k-','linewidth',.5); hold on
 plot([0 0],ax([3 4]),'k-','linewidth',.5); 
 set(gca,'fontsize',20)
 set(gca,'xtick',[-4:2:12])
 title('2-step Backward Difference Formula', 'fontsize',20)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3-step Backward Difference Formula: stability is exterior of boundary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 alpha = [11 -18 9 -2];
 beta  = [6 0 0 0];
 ax = [-4 12 -8 8];
 figure(9), clf
 axis equal, axis(ax)
 fill(ax([1 2 2 1]),ax([3 3 4 4]),fillcolor); hold on
 stabcurve = polyval(alpha,T)./polyval(beta,T);
 obj = fill(real(stabcurve),imag(stabcurve),[1 1 1]);
 axis equal, axis(ax);
 set(obj,'edgecolor',fillcolor)
 plot(ax([1 2]),[0 0],'k-','linewidth',.5); hold on
 plot([0 0],ax([3 4]),'k-','linewidth',.5); 
 set(gca,'fontsize',20)
 set(gca,'xtick',[-4:2:12])
 title('3-step Backward Difference Formula', 'fontsize',20)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4-step Backward Difference Formula: stability is exterior of boundary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 alpha = [25 -48 36 -16 3];
 beta  = [12 0 0 0 0];
 ax = [-4 12 -8 8];
 figure(10), clf
 axis equal, axis(ax)
 fill(ax([1 2 2 1]),ax([3 3 4 4]),fillcolor); hold on
 stabcurve = polyval(alpha,T)./polyval(beta,T);
 obj = fill(real(stabcurve),imag(stabcurve),[1 1 1]);
 axis equal, axis(ax);
 set(obj,'edgecolor',fillcolor)
 plot(ax([1 2]),[0 0],'k-','linewidth',.5); hold on
 plot([0 0],ax([3 4]),'k-','linewidth',.5); 
 set(gca,'fontsize',20)
 set(gca,'xtick',[-4:2:12])
 title('4-step Backward Difference Formula', 'fontsize',20)