% Displays the absolute stability region for a linear multistep method.

% For notation see class notes or Kincaid and Cheney
% for representation of polynomials as vectors in Matlab see help polyval
%
% Fernando Guevara Vasquez 2009

% Euler's method
p = [1 -1]; q =[0 1]; xb=[-2,2]; yb=[-2,2];
% Adams-Bashforth order 3 (three step explicit)
%p = [1, -1, 0, 0]; q=[0,23, -16, 5]/12; xb=[-1,1]; yb=[-1,1];
% Adams-Moulton order 3 (two step implicit)
%p = [1, -1, 0]; q=[5, 8, -1]/12; xb=[-8,5]; yb=[-5,5];


% grid the complex plane
nx=50; ny=50;
x = linspace(xb(1),xb(2),nx); y=linspace(yb(1),yb(2),ny);
[xx,yy]=meshgrid(x,y);
om=xx+1i*yy;
for j=1:ny,
 for i=1:nx,
  maxroot(i,j) = max(abs(roots(p-om(i,j)*q)));
 end;
end;

% draw the contour plot of the magnitude of the largest root of p-omega*q
% for omega in the complex plane
clf;
contour(xx,yy,maxroot,[1 1],'k');
% draw real and imaginary axis
hold on; plot(xb,[0,0],'k'); plot([0,0],yb,'k'); hold off;
axis equal; axis([xb,yb]);
xlabel('real'); ylabel('imaginary');