%StabregionsImplicit
clear;clf
% grid in complex plane
xrm=-2;xrM=14;h=0.05;
yrm=-8;yrM=8;
R = xrm:h:xrM;nx=length(R);
I = yrm:h:yrM;ny=length(I);
% corresponding z values on grid
[x,y]=meshgrid(R,I);
z=x'+1i*y';
for ir=1:nx
    for ii=1:ny
        zz=z(ir,ii);
        if zz~=1,G1(ir,ii)=max(abs(roots([(1-zz) -1])));end                 % BDF1 = BE
        G2(ir,ii)=max(abs(roots([(1-2*zz/3) -4/3 1/3])));                   % BDF2
        G3(ir,ii)=max(abs(roots([(1-6*zz/11) -18/11 9/11 -2/11])));         % BDF3
        G4(ir,ii)=max(abs(roots([(1-12*zz/25) -48/25 36/25 -16/25 3/25]))); % BDF4
   end
end
axis square
contour(R,I,G1',[1 1],'r')
contour(R,I,G2',[1 1],'r')
hold on
contour(R,I,G3',[1 1],'r')
contour(R,I,G4',[1 1],'r')
grid
xlabel('Re(\lambda h)')
ylabel('Im(\lambda h)')
title(' Stability regions for the Backward Difference Schemes ')
axis square


