--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/toolbox/pscont.m	Thu May 07 18:36:24 2015 +0200
@@ -0,0 +1,104 @@
+function [x, y, z, m] = pscont(A, k, npts, ax, levels)
+%PSCONT   Contours and colour pictures of pseudospectra.
+%         PSCONT(A, K, NPTS, AX, LEVELS) plots LOG10(1/NORM(R(z))),
+%         where R(z) = INV(z*I-A) is the resolvent of the square matrix A,
+%         over an NPTS-by-NPTS grid.
+%         NPTS defaults to a SIZE(A)-dependent value.
+%         The limits are AX(1) and AX(2) on the x-axis and
+%                        AX(3) and AX(4) on the y-axis.
+%         If AX is omitted, suitable limits are guessed based on the
+%         eigenvalues of A.
+%         The eigenvalues of A are plotted as crosses `x'.
+%         K determines the type of plot:
+%             K = 0 (default) PCOLOR and CONTOUR
+%             K = 1           PCOLOR only
+%             K = 2           SURFC (SURF and CONTOUR)
+%             K = 3           SURF only
+%             K = 4           CONTOUR only
+%         The contours levels are specified by the vector LEVELS, which
+%         defaults to -10:-1 (recall we are plotting log10 of the data).
+%         Thus, by default, the contour lines trace out the boundaries of
+%         the epsilon pseudospectra for epsilon = 1e-10, ..., 1e-1.
+%         [X, Y, Z, NPTS] = PSCONT(A, ...) returns the plot data X, Y, Z
+%         and the value of NPTS used.
+%
+%         After calling this function you may want to change the
+%         color map (e.g., type COLORMAP HOT - see HELP COLOR) and the
+%         shading (e.g., type SHADING INTERP - see HELP INTERP).
+%         For an explanation of the term `pseudospectra' see PS.M.
+%         When A is real and the grid is symmetric about the x-axis, this
+%         routine exploits symmetry to halve the computational work.
+
+%         Colour pseduospectral pictures of this type are referred to as
+%         `spectral portraits' by Godunov, Kostin, and colleagues.
+%         References:
+%         V. I. Kostin, Linear algebra algorithms with guaranteed accuracy,
+%            Technical Report TR/PA/93/05, CERFACS, Toulouse, France, 1993.
+%         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.
+
+Areal = ~norm(imag(A),1);
+
+if nargin < 5, levels = -10:-1; end
+e = eig(A);
+if nargin < 4
+   ax = cpltaxes(e);
+   if Areal, ax(3) = -ax(4); end  % Make sure region symmetric about x-axis.
+end
+n = max(size(A));
+if nargin < 3, npts = round( min(max(5, sqrt(20^2*10^3/n^3) ), 30)); end
+if nargin < 2, k = 0; end
+
+nptsx = npts; nptsy = npts;
+Ysymmetry = (Areal & ax(3) == -ax(4));
+
+x = seqa(ax(1), ax(2), npts);
+y = seqa(ax(3), ax(4), npts);
+if Ysymmetry                    % Exploit symmetry about x-axis.
+   nptsy = ceil(npts/2);
+   y1 = y;
+   y = y(1:nptsy);
+end
+
+[xx, yy] = meshgrid(x,y);
+z = xx + sqrt(-1)*yy;
+I = eye(n);
+Smin = zeros(nptsy, nptsx);
+
+for j=1:nptsx
+    for i=1:nptsy
+        Smin(i,j) = min( svd( z(i,j)*I-A ) );
+    end
+end
+
+z = log10( Smin + eps );
+if Ysymmetry
+   z = [z; z(nptsy-rem(npts,2):-1:1,:)];
+   y = y1;
+end
+
+if k == 0 | k == 1
+   pcolor(x, y, z); hold on
+elseif k == 2
+   surfc(x, y, z); hold on
+elseif k == 3
+   surf(x, y, z); hold on
+end
+
+if k == 0 | k == 4
+   contour(x, y, z, levels); hold on
+end
+
+if k ~= 2 & k ~= 3
+   if k == 0 | k == 1
+      s = 'k';   % Black.
+   else
+      s = 'w';   % White.
+   end
+   plot(real(e),imag(e),['x' s]);
+end
+
+axis('square');
+hold off