% plot real pseudo spectrum
function []=rpseudo(a,e,m)
% inputs:
%   a = matrix
%   e = size of perturbation
%   m = number of perturbed versions
[n,n]=size(a);
d=0;
d(1:n,1) = eig(a);
for i=1:m,
  [q,r]=qr(randn(n,n));
  ap=q'*(a+ e*(2*rand(n,n)-1))*q;
  d(1:n,i+1)=eig(ap);
end
scl = max(max(abs(d)));
sclr = max(max(abs(real(d))))*1.2 + 1e-15*scl;
scli = max(max(abs(imag(d))))*1.2 + 1e-15*scl;
figure(1),
hold off, clg
plot(real(d(:,1)),imag(d(:,1)),'*r')
title('Real Pseudo-spectrum')
axis([-sclr sclr -scli scli]),
hold on
for i=1:m
  plot(real(d(:,i+1)),imag(d(:,i+1)),'+w')
end