
th=0:.001:2*pi; th=th';
ii=sqrt(-1); ith=ii*th; ei=exp(ith); ep=1.e-13;

E = [ ei 1+0*ei 1./ei 1./(ei.*ei) 1./(ei.*ei.*ei)];


ab0   = [1 0.0 0.0 0. 0.]';
ab1   = [0 1.0 0.0 0. 0.]';
ab2   = [0 1.5 -.5 0. 0.]';
ab3   = [0 23./12. -16./12. 5./12. 0.]';

du    = [1. -1. 0. 0. 0.]';

ymax  = 2; yaxis = [-ymax*ii ymax*ii]'; xaxis = [-3.0+ep*ii 3.0+ep*ii]';

hold off; 
plot (yaxis,'k-'); axis square; axis([-ymax-.5 ymax-.5 -ymax ymax]); hold on;
plot (xaxis,'k-'); 

ldteb     = (E*du)./(E*ab0);         plot (ldteb  ,'g-');
ldtef     = (E*du)./(E*ab1);         plot (ldtef  ,'b-');
ldtab3    = (E*du)./(E*ab3);         plot (ldtab3 ,'r-');

title('Stability Regions for EF, EB, AB3')
xlabel('Re \lambda\Delta t')
ylabel('Im \lambda\Delta t')

