function y = ps(A, m, tol, rl, marksize)
%PS     Dot plot of a pseudospectrum.
%       PS(A, M, TOL, RL) plots an approximation to a pseudospectrum
%       of the square matrix A, using M random perturbations of size TOL.
%       M defaults to a SIZE(A)-dependent value and TOL to 1E-3.
%       RL defines the type of perturbation:
%         RL =  0 (default): absolute complex perturbations of 2-norm TOL.
%         RL =  1:           absolute real perturbations of 2-norm TOL.
%         RL = -1:           componentwise real perturbations of size TOL.
%       The eigenvalues of A are plotted as crosses `x'.
%       PS(A, M, TOL, RL, MARKSIZE) uses the specified marker size instead
%       of a size that depends on the figure size, the matrix order, and M.
%       If MARKSIZE < 0, the plot is suppressed and the plot data is returned
%       as an output argument.
%       PS(A, 0) plots just the eigenvalues of A.

%       For a given TOL, the pseudospectrum of A is the set of
%       pseudo-eigenvalues of A, that is, the set
%       { e : e is an eigenvalue of A+E, for some E with NORM(E) <= TOL }.
%
%       Reference:
%       L.N. Trefethen, Pseudospectra of matrices, in D.F. Griffiths and
%            G.A. Watson, eds, Numerical Analysis 1991, Proceedings of the 14th
%            Dundee Conference, vol. 260, Pitman Research Notes in Mathematics,
%            Longman Scientific and Technical, Essex, UK, 1992, pp. 234-266.

if diff(size(A)), error('Matrix must be square.'), end
n = max(size(A));

if nargin < 5, marksize = 0; end
if nargin < 4, rl = 0; end
if nargin < 3, tol = 1e-3; end
if nargin < 2, m = max(1, round( 25*exp(-0.047*n) )); end

if m == 0
   e = eig(A);
   ax = cpltaxes(e);
   plot(real(e), imag(e), 'xg')
   axis(ax);
   axis('square');
   return
end

x = zeros(m*n,1);
i = sqrt(-1);

for j = 1:m
   if rl == -1     % Componentwise.
      dA = -ones(n) + 2*rand(n);   % Uniform random numbers on [-1,1].
      dA = tol * A .* dA;
   else
      if rl == 0   % Complex absolute.
         dA = randn(n) + i*randn(n);
      else         % Real absolute.
         dA = randn(n);
      end
      dA = tol/norm(dA)*dA;
   end
   e = eig(A + dA);
   x((j-1)*n+1:j*n) = e;
end

if marksize >= 0

   ax = cpltaxes(x);
   h = plot(real(x),imag(x),'.');
   axis(ax);
   axis('square');

   % Next block adapted from SPY.M.
   if marksize == 0
      units = get(gca,'units');
      set(gca,'units','points');
      pos = get(gca,'position');
      nps = 2.4*sqrt(n*m);  % Factor based on number of pseudo-ei'vals plotted.
      myguess = round(3*min(pos(3:4))/nps);
%      [nps myguess]
      marksize = max(1,myguess);
      set(gca,'units',units);
   end
   set(h,'markersize',marksize);
%   set(h,'linemarkersize',marksize);

   hold on
   e = eig(A);
   plot(real(e),imag(e),'xw');
   set(h,'markersize',marksize);
   hold off

else

  y = x;

end