%
%  general script for stability diagram
%

clear

%  define complex grid

qr = -2:0.05:2;
qi = -2:0.05:2;

%  loop for computing amplification factors

for Nr = 1:size(qr,2)
for Ni = 1:size(qi,2)

  qh0 = qr(Nr) + i*qi(Ni);
  qh(Ni,Nr) = qh0;

  lambda(Ni,Nr) = abs(1 + qh0 + qh0^2/2 + qh0^3/2);

end
end

%  make plot

figure(1);  clf
contour(qr,qi,lambda,'b')
hold on
contour(qr,qi,lambda,[1 1],'r')
plot([0 0],[-10 10],'k--');  plot([-10 10],[0 0],'k--')
title('\bf stability diagram')
xlabel('\bf real axis');  ylabel('\bf imag axis')

figure(2);  clf
mid = (size(qi,2)+1)/2;
midlam = lambda(mid,:);

plot(qr,midlam,'rx')
hold on
plot(qr,midlam./(midlam>1),'bx')
title('\bf zeros for real q')
xlabel('\bf q\Delta t axis');  ylabel('\bf |amplification|')