function ezgraph3(surfstyle,varargin)
%EZGRAPH3 Helper function for easy surface plotters.
%   EZGRAPH3(PLOTFUN,f) plots z=f(x,y) using the plotting function PLOTFUN.
%   PLOTFUN is any function like SURF that takes (x,y,z) as inputs.
%
%   EZGRAPH3(PLOTFUN,f,[xmin xmax ymin ymax]) plots over the specified domain.
%   EZGRAPH3(PLOTFUN,...,DOMSTYLE) DOMSTYLE can be 'rect' or 'circ'.
%   EZGRAPH3(PLOTFUN,...,Npts) uses a grid of Npts-by-Npts points.
%
%   EZGRAPH3(PLOTFUN,x,y,z) plots the paramteric surface x=x(s,t), y=y(s,t),
%   and z=z(s,t) over the square -2*pi < s < 2*pi and -2*pi < t < 2*pi.
%
%   EZGRAPH3(PLOTFUN,x,y,z,[smin,smax,tmin,tmax]) or EZGRAPH3(x,y,z,[a,b])
%   uses the specified domain. 
%
%   See also EZPLOT, EZPLOT3, EZPOLAR, EZMESH, EZMESHC, EZSURF,
%            EZSURFC, EZCONTOUR. EZCONTOURF.

%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.2 $  $Date: 1998/07/10 16:41:13 $
%   P. J. Costa and Clay Thompson

% Set defaults
NA = nargin - 1;
rotate3d on;
fig = gcf;

[ezflag,f,domain,domstyle,Npts] = parse_inputs(varargin{:});

if ezflag % Use ezplot3s
   ezplot3s(Npts,f{:},domain,surfstyle);
elseif ~ezflag & isequal(domstyle,'rect') 
   % Plot a surface z = f(x,y) over a rectangle [x0,xf] x [y0,yf].
   surfplot(f,domain,surfstyle,fig,Npts,NA);
elseif isequal(domstyle,'circ')
   % Plot a surface z = f(x,y) over a circle by letting x = r*cos(t),
   % y = r*sin(t), z = f(r*cos(t),r*sin(t)), r in [0,rm] and t in [-pi,pi].
   a = domain(1); b = domain(2); c = domain(3); d = domain(4);
   Var = symvar(f);
   if isempty(Var); Var = {'x','y'}; end
   Cen = [a+b,c+d]/2;
   % Change variables to (r,t)
   x = ['(r*cos(t) + ' num2str(Cen(1)) ')'];
   y = ['(r*sin(t) + ' num2str(Cen(1)) ')'];
   if isa(f,'inline')
     z = inline(idrep(idrep(char(f),Var{1},x),Var{2},y),'r','t');
   else
     z = idrep(idrep(char(f),Var{1},x),Var{2},y);
   end
   rm = sqrt((Cen(1) - b)^2 + (Cen(2) - d)^2);
   domain = [0,rm,-pi,pi];
   ezplot3s(76,x,y,z,domain,surfstyle);
   title(texlabel(f));
end 

%--------------------------------------------
function c = commas(s)
%COMMAS  COMMAS(s) is true for commas not inside parentheses.
p = cumsum((s == '(') - (s == ')'));
c = (s == ',') & (p == 0);

%-----------------------------------------
function R = dom(v)
%DOM creates a 1-by-4 class double vector whose elements
%   correspond to the endpoints of a rectangle.  If the
%   length of the input vector v is 1, then the rectangle R
%   is [-|v|,|v|,-|v|,-|v|].  If the length of the input
%   is 2, then R = [-v(1),v(2),-v(1),v(2)].  If the length
%   of v is 3, the R = [v(1),v(2),-|v(3)|,|v(3)|] and if the
%   length of of v is 4, then R = v.  Otherwise,
%   R = [-|v(1)|,|v(1)|,-|v(1)|,-|v(1)|]
switch length(v)
   case 0
      R = [-2*pi, 2*pi, -2*pi, 2*pi];
   case 1
      R = double([-abs(v),abs(v),-abs(v),abs(v)]);
   case 2
      R = double([v(1),v(2),v(1),v(2)]);
   case 3
      R = double([-v(1),v(2),-abs(v(3)),abs(v(3))]);
   case 4
      R = double(v);
   otherwise
      R = double([-abs(v(1)), abs(v(1)), -abs(v(1)), abs(v(1))]);
end

%-------------------------

function ezplot3s(Npts,varargin)
% EZPLOT3S  Easy to use 3-d parametric curve or surface plotter.
%   EZPLOT3S(Npts,x,y,z) plots the surfaces S(t,s) = (x(t,s),y(t,s),z(t,s))
%   using a grid of Npts-by-Npts points.  The default region of surfaces is
%   the square [-2*pi,2*pi] X [-2*pi,2*pi].
%
%   EZPLOT3S(Npts,x,y,z,[smin,smax],[tmin,tmax]), plots the surface
%   S(t,s) = (x(t,s),y(t,s),z(t,s)) for smin < s < smax, tmin < t < tmax.
%   In the case thatðx,y,z are not functions of s and t 
%   (rather, say u and v), then the domain endpoints [umin,umax] 
%   [vmin,vmax] are given alphabetically.
%
%   EZPLOT3S(Npts,x,y,z,[a,b]) plots S(t,s) for a < s < b, a < t < b.
%   
%   Examples
%     syms t s
%     ezplot3s(50,s*cos(t),s*sin(t),t)
%     ezplot3s(50,s*cos(t),s*sin(t),s)
%     ezplot3s(50,exp(-s)*cos(t),exp(-s)*sin(t),t,[0,8],[0,2*pi])
%     ezplot3s(50,cos(s)*cos(t),cos(s)*sin(t),sin(s),[0,pi/2],[0,3*pi/2])

% At the beginning of the calculation, set the pointer to an hour glass.
set(gcf,'Pointer','watch');
NA = length(varargin);
x = varargin{1}; y = varargin{2}; z = varargin{3};

% Determine the vairables in that define the parametrization
varsx = symvar(x);
varsy = symvar(y);
varsz = symvar(z);
svars = union(union(varsx,varsy),varsz);
% Determine the number of variables in r = (x,y,z)
N = length(svars);
if isequal(svars,'{''}')
   svars = {};
end
% Determine the surface style
if isa(varargin{end},'char')
   surfstyle = varargin{end};
   L = length(varargin)-1;
else
   surfstyle = 'surf';
   L = length(varargin);
end

% Determine the domains in t and s:
switch L
   case 3
      T = linspace(-2*pi,2*pi,Npts);
      S = linspace(-2*pi,2*pi,Npts);
   case 4
      if ( prod(size(varargin{4})) == 4 )
         T = linspace(varargin{4}(3),varargin{4}(4),Npts);
         S = linspace(varargin{4}(1),varargin{4}(2),Npts);
      else
         T = linspace(varargin{4}(1),varargin{4}(2),Npts);
         S = T;
      end
   case 5
      T = linspace(varargin{5}(1),varargin{5}(2),Npts);
      S = linspace(varargin{4}(1),varargin{4}(2),Npts);
end

switch N
   case 0 % ezsurf('2','3','4')
      plot3(eval(x),eval(y),eval(z),'b.','MarkerSize',24);
   case 1 
      error('To plot parametric curves, use EZPLOT3.');
   case 2 % Parametrized surface
      V = svars;
      if isempty(V), V{1} = 't'; V{2} = 's'; end
      [S,T] = meshgrid(S,T);
      X = eval(vectorize(idrep(idrep(char(x),V{1},'S'),V{2},'T')));
      % Special processing in case X, Y, or Z are 1-by-1 syms/doubles.
      if prod(size(X)) == 1
         Xmin = 0.9*X; Xmax = 1.1*X;
         X = X*ones(size(S));
      else
         Xmin = min(min(X)); Xmax = max(max(X));
      end
      Y = eval(vectorize(idrep(idrep(char(y),V{1},'S'),V{2},'T')));
      if prod(size(Y)) == 1
         Ymin = 0.9*Y; Ymax = 1.1*Y;
         Y = Y*ones(size(S));
      else
         Ymin = min(min(Y)); Ymax = max(max(Y));
      end
      Z = eval(vectorize(idrep(idrep(char(z),V{1},'S'),V{2},'T')));
      if prod(size(Z)) == 1
         Zmin = 0.9*Z; Zmax = 1.1*Z;
         Z = Z*ones(size(S));
      else
         Zmin = min(min(Z)); Zmax = max(max(Z));
      end

      [S,T,X,xmin,xmax] = gridfilt(S,T,X,Npts,NA,x);
      [S,T,Y,ymin,ymax] = gridfilt(S,T,Y,Npts,NA,y);
      [S,T,Z,zmin,zmax] = gridfilt(S,T,Z,Npts,NA,z);

      % Plot the paramterically defined surface.
      feval(surfstyle,X,Y,Z);
      xlabel('x'); ylabel('y'); zlabel('z');
      % Set the appropriate limits in x, y, and z.
      Xmin = 1.2*min(Xmin,xmin); Xmax = 1.2*max(Xmax,xmax);
      Ymin = 1.2*min(Ymin,ymin); Ymax = 1.2*max(Ymax,ymax);
      % Change Zmin and Zmax to get entire plot within the axis frame.
      Zmin = 1.2*min(Zmin,zmin); Zmax = 1.2*max(Zmax,zmax);
      if isequal(surfstyle,'contour') | isequal(surfstyle,'contourf')
         axis([Xmin,Xmax,Ymin,Ymax]);
      else
         axis([Xmin,Xmax,Ymin,Ymax,Zmin,Zmax]);
      end
      title(['x = ' texlabel(x), ', y = ' texlabel(y), ', z = ' texlabel(z)]);
   otherwise
      error('Cannot plot in 4 dimensions')
end

% At the end of the calculation, reset the pointer to an arrow.
set(gcf,'Pointer','arrow');

%----------------------------------------

function [X,Y,u,umin,umax] = gridfilt(X,Y,u,N,NA,F)
%GRIDFILT filters the computational grid to remove singularities and
%   imaginary (complex) numbers.
%   GRIDFILT(X,Y,Z,N,NA) returns a surface grid [X,Y], u = F(X,Y) whose
%   singularities and complex number values have been set to NaNs.  Here
%   X,Y, and u represent the points in a surface, N = number of points
%   in the grid, NA = nargin, and F is the inline function that defines
%   u:  u = F(X,Y).

t = (0:N-1)/(N-1);
% Reduce to an "interesting" domain.
if (N > 1) & (NA <= 3)
   dx = X(2)-X(1);
   if (dx == 0), dx = 0.025; end
   dy = Y(2)-Y(1);
   if (dy == 0), dy = 0.025; end
   if prod(size(u)) == 1
      dudx = 0; dudy = 0;
   else
      [dudx,dudy] = gradient(u,dx,dy);
   end
   
   k = find(abs(dudx)+abs(dudy) > .01);
   if isempty(k)
      xmin = min(X(:)); xmax = max(X(:));
      ymin = min(Y(:)); ymax = max(Y(:));
   else
      xmin = min(X(k)); xmax = max(X(k));
      ymin = min(Y(k)); ymax = max(Y(k));
   end
   if xmin < floor(4*xmin)/4 + dx, xmin = floor(4*xmin)/4; end
   if xmax > ceil(4*xmax)/4 - dx, xmax = ceil(4*xmax)/4; end
   X = xmin + t*(xmax-xmin);
   if ymin < floor(4*ymin)/4 + dy, ymin = floor(4*ymin)/4; end
   if ymax > ceil(4*ymax)/4 - dy, ymax = ceil(4*ymax)/4; end
   Y = ymin + t*(ymax-ymin);
   [X,Y] = meshgrid(X,Y);
   u = feval(F,X,Y);
   Umin = min(u(:)); Umax = max(u(:));
   k = find(abs(imag(u)) > 1.e-6*abs(real(u)));
   if any(k), u(k) = NaN; end
end

% Determine u scale so that "most" of the u values
% are in range, but singularities are off scale.
u = real(u);
uu = sort(u(isfinite(u)));
N = length(uu);
if N > 16
   del = uu(fix(15*N/16)) - uu(fix(N/16));
   umin = max(uu(1)-del/16,uu(fix(N/16))-del);
   umax = min(uu(N)+del/16,uu(fix(15*N/16))+del);
elseif N > 0
   umin = uu(1);
   umax = uu(N);
else
   umin = 0;
   umax = 0;
end
if umin == umax, umin = umin-1; umax = umax+1; end

% Eliminate vertical lines at discontinuities.
[p,q] = find(abs(del2(u)) > .4);
U = u;
if (~isempty(p) & ~isempty(q)), U(p,q) = NaN; end
unan = isnan(U) | u > umax | u < umin;
%  If everything gets set to a NaN proceed via the ezplot/Moler
%  filtering . . .
if ( sum(sum(unan))/prod(size(U)) > 0.50 )
   [nr,nc] = size(u);
   % First, search along the rows . . .
   for j = 1:nr
     k = 2:nc;
     krow = find( (u(j,k) > umax/3 & u(j,k-1) < umin/3) | ...
                 (u(j,k) < umin/3 & u(j,k-1) > umax/3) );
     if any(krow), u(j,krow) = NaN; end
   end
   % . . . then search along the columns.
   for j = 1:nc
     k = 2:nr;
     kcol = find( (u(k,j) > umax/3 & u(k-1,j) < umin/3) | ...
                 (u(k,j) < umin/3.4 & u(k-1,j) > umax/3) );
     if any(kcol), u(kcol,j) = NaN; end
   end
% . . . otherwise use the Laplacian = 0 filtering
elseif (~isempty(p) & ~isempty(q))
   u(p,q) = NaN;
else
   u = u;
end

% Limit data to umin and umax to "fake" clipping
if (~isempty(p) & ~isempty(q))
   k = find(u > umax | u < umin);
   u(k) = NaN;
end

%-------------------------

function s = idrep(s1,s2,s3)
%IDREP Replaces all occurences of s2 in the string s1 by s3.
%   For example, idrep(x*exp(x*y),x,W) returns W*exp(W*y).

% Take care of trivial cases:
%   s2 > s1
%   s1 is empty
%   s2 is empty
if (length(s2) > length(s1)) | isempty(s2) | isempty(s1),
  s=s1;
  return;
end;

% Place non-identifiers at the beginning and end of the string.
s=[]; s1 = [ '#' s1 '#' ];
% Find all occurrences of s2 in s1
s2len=length(s2);
s2pos=findstr(s1,s2);

% Build resulting string
s1pos=1;
for i=1:length(s2pos),
  if ~isalphanum(s1(s2pos(i)-1)) & ~isalphanum(s1(s2pos(i)+s2len))
     s=[s,s1(s1pos:s2pos(i)-1),s3];
     s1pos=s2pos(i)+s2len;
  end
end;
s=[s,s1(s1pos:length(s1))];
% Remove the #'s at the beginning and end of s.
s(1) = []; s(end) = [];

%-------------------------

function B = isalphanum(S)
%ISALPHANUM True for alpha-numeric characters.
%   For a string S, ISALPHANUM returns 1 for alpha-numeric
%   characters (or underscores) and 0 otherwise.
%   Example:
%      S = 'x*exp(x - y) + cosh(x*z^2)';
%      isalphanum(S)  returns
%      (1,0,1,1,1,0,1,0,0,0,1,0,0,0,0,1,1,1,1,0,1,0,1,0,1,0)

B = isletter(S) | (S >= '0' & S <= '9') | (S == '_');

%-----------------------------

function tf = isfun(f)
% ISFUN Determines whether an expression is a character
%   string or on inline function.

tf = isa(f,'char') | isa(f,'inline');

%--------------------------------------------------

function [ezflag,f,domain,domstyle,Npts] = parse_inputs(varargin)
% PARSE_INPUTS Parses the inputs to EZGRAPH3.

ezflag = 0;  % EZFLAG = 0 when ezplot3s is NOT in use.
domain = [-2*pi, 2*pi, -2*pi, 2*pi];
Npts = 60;
domstyle = 'rect';
f = varargin{1};
V = symvar(f);
if isempty(V), V = {'x','y'}; end
% If the function is an M-file (i.e., 'peaks'), assign
% the variables {'x','y'}.
if exist(char(f)) == 2 & length(V) == 1
   V = {'x','y'};
end
% If the function is an M-file, make it an inline of the form
% foo(x,y).
if exist(char(f)) == 2
   f = inline([f '(' V{1} ',' V{2} ')'],V{:});
end

switch length(varargin)
case 0
   error('Mathematical expression required.');
case 1
   % ezgraph3(surfstyle,f) - use the default values
case 2
   % ezgraph3(surfstyle,f,domain)
   if isfun(varargin{1}) & isa(varargin{2},'double') & ...
      (length(varargin{2}) > 1)
      domain = dom(varargin{2});

   % ezgraph3(surfstyle,f,Npts)
   elseif isfun(varargin{1}) & isa(varargin{2},'double') & ...
      length(varargin{2}) == 1
      Npts = varargin{2};

   % ezgraph3(surfstyle,f,domstyle)
   elseif isfun(varargin{1}) & isa(varargin{2},'char')
      domstyle = varargin{2};
   end
case 3
   % ezgraph3(surfstyle,x,y,z)
   if isfun(varargin{1}) & isfun(varargin{2}) & isfun(varargin{3})
      ezflag = 1;
      x = varargin{1}; y = varargin{2}; z = varargin{3};
      Npts = 50;

   % ezgraph3(surfstyle,f,domain(x,y),domstyle)
   elseif isfun(varargin{1}) & isa(varargin{2},'double') & ...
          (length(varargin{2}) > 1) & isa(varargin{3},'char')
      domain = dom(varargin{2});
      domstyle = varargin{3};

   % ezgraph3(surfstyle,f,Npts,domstyle)
   elseif isfun(varargin{1}) & isa(varargin{2},'double') & ...
          (length(varargin{2}) == 1) & isa(varargin{3},'char')
      Npts = varargin{2};
      domstyle = varargin{3};

   % ezgraph3(surfstyle,f,domain(x),domain(y))
   elseif isfun(varargin{1}) & isa(varargin{2},'double') & ...
          isa(varargin{3},'double') & ...
          (length(varargin{2}) > 1) & (length(varargin{3}) > 1)
      domain = dom([varargin{2}, varargin{3}]);

   % ezgraph3(surfstyle,f,domain(x,y),Npts)
   elseif isfun(varargin{1}) & isa(varargin{2},'double') & ...
          isa(varargin{3},'double') & ...
          (length(varargin{2}) > 1) & (length(varargin{3}) == 1)
      domain = dom(varargin{2});
      Npts = varargin{3};
   else
      error('Incorrect inputs.');
   end
case 4
   % ezgraph3(surfstyle,x,y,z,domain(s,t))
   if isfun(varargin{1}) & isfun(varargin{2}) & ...
      isfun(varargin{3}) & ...
      isa(varargin{4},'double') & (length(varargin{4}) > 1)
      ezflag = 1;
      x = varargin{1}; y = varargin{2}; z = varargin{3};
      domain = dom(varargin{4});
      Npts = 50;

   % ezgraph3(surfstyle,x,y,z,Npts)
   elseif isfun(varargin{1}) & isfun(varargin{2}) & ...
          isfun(varargin{3}) & ...
          isa(varargin{4},'double') & (length(varargin{4}) == 1)
      ezflag = 1;
      x = varargin{1}; y = varargin{2}; z = varargin{3};
      Npts = varargin{4};

   % ezgraph3(surfstyle,x,y,z,domstyle)
   elseif isfun(varargin{1}) & isfun(varargin{2}) & ...
          isfun(varargin{3}) & isa(varargin{4},'char')
      ezflag = 1;
      x = varargin{1}; y = varargin{2}; z = varargin{3};
      domstyle = varargin{4};
      Npts = 50;

   % ezgraph3(surfstyle,f,domain(x),domain(y),domstyle)
   elseif isfun(varargin{1}) & isa(varargin{4},'char') & ...
          isa(varargin{2},'double') & (length(varargin{2}) > 1) & ...
          isa(varargin{3},'double') & (length(varargin{3}) > 1)
      domain = dom([varargin{2},varargin{3}]);
      domstyle = varargin{4};

   % ezgraph3(surfstyle,f,domain(x),domain(y),Npts)
   elseif isfun(varargin{1}) & ...
          isa(varargin{2},'double') & (length(varargin{2}) > 1) & ...
          isa(varargin{3},'double') & (length(varargin{3}) > 1) & ...
          isa(varargin{4},'double') & (length(varargin{4}) == 1)
      domain = dom([varargin{2},varargin{3}]);
      Npts = varargin{4};

   % ezgraph3(surfstyle,f,domstyle,domain(x,y),Npts)
   elseif isfun(varargin{1}) & isa(varargin{2},'char') & ...
          isa(varargin{3},'double') & (length(varargin{3}) > 1) & ...
          isa(varargin{4},'double') & (length(varargin{4}) == 1)
      domstyle = varargin{2};
      domain = dom(varargin{3});
      Npts = varargin{4};

   % ezgraph3(surfstyle,f,domain(x,y),Npts,domstyle)
   elseif isfun(varargin{1}) & isa(varargin{2},'double') & ...
          (length(varargin{2}) > 1) & isa(varargin{4},'char') & ...
          isa(varargin{3},'double') & (length(varargin{3}) == 1)
      domain = dom(varargin{2});
      Npts = varargin{3};
      domstyle = varargin{4};
   else
      error('Incorrect inputs.');
   end
case 5
   % ezgraph3(surfstyle,x,y,z,domain(s,t),Npts)
   if isfun(varargin{1}) & (length(varargin{4}) > 1) & ...
         (length(varargin{5}) == 1)
      ezflag = 1;
      x = varargin{1}; y = varargin{2}; z = varargin{3};
      domain = dom(varargin{4});
      Npts = varargin{5};

   % ezgraph3(surfstyle,x,y,z,domain(s),domain(t))
   elseif isfun(varargin{1}) & (length(varargin{4}) == 2) & ...
          (length(varargin{5}) == 2)
      ezflag = 1;
       x = varargin{1}; y = varargin{2}; z = varargin{3};
      domain = dom([varargin{4},varargin{5}]);
      Npts = 50;
   else
      error('Incorrect inputs.');
   end
case 6
   % ezgraph3(surfstyle,x,y,z,domain(s),domain(t),Npts)
   if isfun(varargin{1}) & (length(varargin{4}) == 2) & ...
          (length(varargin{5}) == 2)
      ezflag = 1;
      x = varargin{1}; y = varargin{2}; z = varargin{3};
      domain = dom([varargin{4},varargin{5}]);
      Npts = varargin{6};
   end
otherwise
   f = varargin{1};
   domain = dom(varargin{2});
   domstyle = varargin{3};
   Npts = varargin{4};
end

if ezflag, % Package x,y,z into f.
  f = {x,y,z};
end

%-------------------------

function zlim = surfplot(f,domain,surfstyle,fig,N,NA)
%SURFPLOT Plots the surface z = f associated with the symbolic expression f
%  of two symbolic variables over the DOMAIN with a surface style SURFSTYLE,
%  figure number FIG, and N-by-N grid. NA is the number of input arguments 
%  from ezgraph3.
%  SURFPLOT(x^2 + y^2,[0,1,2,3],surf,3,50,1) plots the surface z = x^2 + y^2,
%  for 0 < x < 1 and 2 < y < 3 via surf in figure(3) over a 50-by-50 grid.

% At the beginning of the calculation, set the pointer to an hour glass.
set(gcf,'Pointer','watch');

F = vectorize(fcnchk(f));

if ischar(F) & exist(F)==2 & nargin(F)==2
  var = {'x','y'};
else
  var = argnames(F);
end

% Make sure we have two variables
switch(length(var))
case 0
   var = {'x','y'};
case 1
  if isequal(var,'y')
    var = {'x', var{1}};
  else
    var = {var{1}, 'y'};
  end
  F = inline(char(F),var{:});
case 2
  var = var;
otherwise
  error(['The expression ' char(f) ' must only have 2 symbolic variables']);
end
    
% Sample on the initial rectangle with X = xmin + t*(xmax-xmin),
% Y = ymin + t*(ymax-ymin) and t = [0, 1/N, 2/N, ... , (N-1)/N];
t = (0:N-1)/(N-1);
xmin = min(domain(1),domain(2));
xmax = max(domain(1),domain(2));
X = xmin + t*(xmax - xmin);
ymin = min(domain(3),domain(4));
ymax = max(domain(3),domain(4));
Y = ymin + t*(ymax - ymin);
[X,Y] = meshgrid(X,Y);
u = feval(F,X,Y);
Umin = min(u(:)); Umax = max(u(:));
k = find(abs(imag(u)) > 1.e-6*abs(real(u)));
if any(k), u(k) = NaN; X(k) = NaN; Y(k) = NaN; N = length(Y); end

[X,Y,u,umin,umax] = gridfilt(X,Y,u,N,NA,F);

% Plot the function
if (isnan(xmin) | isnan(xmax) | isnan(ymin) | isnan(ymax) )
   error([char(f) ' must be real-valued only.'])
end
figure(fig);
if prod(size(u)) == 1, u = u*ones(size(X)); end
feval(surfstyle,X,Y,u);

% Set the appropriate limits in x, y, and z.
zlim = get(gca,'zlim');
xmin = max(xmin,min(min(X))); xmax = min(xmax,max(max(X)));
ymin = max(ymin,min(min(Y))); ymax = min(ymax,max(max(Y)));

% Change umin and umax to get entire plot within the axis frame.
if abs(Umin - umin) > 0.05*abs(Umin) | Umin == -Inf, umin = zlim(1); end
if abs(Umax - umax) > 0.05*abs(Umax) | Umax == Inf, umax = zlim(2); end
if isequal(surfstyle,'contour') | isequal(surfstyle,'contourf')
   axis([xmin xmax ymin ymax]);
   grid off;
else
   axis([xmin xmax ymin ymax umin umax]);
   grid on;
 end

if ~isnan(var{1}), xlabel(texlabel(var{1})); end
if ~isnan(var{2}), ylabel(texlabel(var{2})); end

% Change the title to TeX format.
title(texlabel(f));

% At the end of the calculation, reset the pointer to an arrow.
set(gcf,'Pointer','arrow');
