--- a/scripts/plot/ezcontour.m	Tue Jun 04 19:41:21 2013 +0200
+++ b/scripts/plot/ezcontour.m	Wed Jul 17 10:09:44 2013 -0700
@@ -20,23 +20,28 @@
 ## @deftypefn  {Function File} {} ezcontour (@var{f})
 ## @deftypefnx {Function File} {} ezcontour (@dots{}, @var{dom})
 ## @deftypefnx {Function File} {} ezcontour (@dots{}, @var{n})
-## @deftypefnx {Function File} {} ezcontour (@var{h}, @dots{})
+## @deftypefnx {Function File} {} ezcontour (@var{hax}, @dots{})
 ## @deftypefnx {Function File} {@var{h} =} ezcontour (@dots{})
 ##
-## Plot the contour lines of a function.  @var{f} is a string, inline function
-## or function handle with two arguments defining the function.  By default the
-## plot is over the domain @code{-2*pi < @var{x} < 2*pi} and @code{-2*pi <
-## @var{y} < 2*pi} with 60 points in each dimension.
+## Plot the contour lines of a function.
+## 
+## @var{f} is a string, inline function, or function handle with two arguments
+## defining the function.  By default the plot is over the meshed domain
+## @code{-2*pi <= @var{x} | @var{y} <= 2*pi} with 60 points in each dimension.
 ##
 ## If @var{dom} is a two element vector, it represents the minimum and maximum
-## value of both @var{x} and @var{y}.  If @var{dom} is a four element vector,
-## then the minimum and maximum value of @var{x} and @var{y} are specify
-## separately.
+## values of both @var{x} and @var{y}.  If @var{dom} is a four element vector,
+## then the minimum and maximum values are @code{[xmin xmax ymin ymax]}.
 ##
 ## @var{n} is a scalar defining the number of points to use in each dimension.
 ##
+## If the first argument is an axis handle, @var{hax}, then plot into this
+## axis rather than the current axis handle returned by @code{gca}.
+##
 ## The optional return value @var{h} is a graphics handle to the created plot.
 ##
+## Example:
+##
 ## @example
 ## @group
 ## f = @@(x,y) sqrt (abs (x .* y)) ./ (1 + x.^2 + y.^2);
@@ -44,19 +49,19 @@
 ## @end group
 ## @end example
 ##
-## @seealso{ezplot, ezcontourf, ezsurfc, ezmeshc}
+## @seealso{contour, ezcontourf, ezplot, ezmeshc, ezsurfc}
 ## @end deftypefn
 
-function retval = ezcontour (varargin)
+function h = ezcontour (varargin)
 
-  [h, needusage] = __ezplot__ ("contour", varargin{:});
+  [htmp, needusage] = __ezplot__ ("contour", varargin{:});
 
   if (needusage)
     print_usage ();
   endif
 
   if (nargout > 0)
-    retval = h;
+    h = htmp;
   endif
 
 endfunction

--- a/scripts/plot/ezcontourf.m	Tue Jun 04 19:41:21 2013 +0200
+++ b/scripts/plot/ezcontourf.m	Wed Jul 17 10:09:44 2013 -0700
@@ -20,23 +20,28 @@
 ## @deftypefn  {Function File} {} ezcontourf (@var{f})
 ## @deftypefnx {Function File} {} ezcontourf (@dots{}, @var{dom})
 ## @deftypefnx {Function File} {} ezcontourf (@dots{}, @var{n})
-## @deftypefnx {Function File} {} ezcontourf (@var{h}, @dots{})
+## @deftypefnx {Function File} {} ezcontourf (@var{hax}, @dots{})
 ## @deftypefnx {Function File} {@var{h} =} ezcontourf (@dots{})
 ##
-## Plot the filled contour lines of a function.  @var{f} is a string, inline
-## function or function handle with two arguments defining the function.  By
-## default the plot is over the domain @code{-2*pi < @var{x} < 2*pi} and
-## @code{-2*pi < @var{y} < 2*pi} with 60 points in each dimension.
+## Plot the filled contour lines of a function.
+## 
+## @var{f} is a string, inline function, or function handle with two arguments
+## defining the function.  By default the plot is over the meshed domain
+## @code{-2*pi <= @var{x} | @var{y} <= 2*pi} with 60 points in each dimension.
 ##
 ## If @var{dom} is a two element vector, it represents the minimum and maximum
-## value of both @var{x} and @var{y}.  If @var{dom} is a four element vector,
-## then the minimum and maximum value of @var{x} and @var{y} are specify
-## separately.
+## values of both @var{x} and @var{y}.  If @var{dom} is a four element vector,
+## then the minimum and maximum values are @code{[xmin xmax ymin ymax]}.
 ##
 ## @var{n} is a scalar defining the number of points to use in each dimension.
 ##
+## If the first argument is an axis handle, @var{hax}, then plot into this
+## axis rather than the current axis handle returned by @code{gca}.
+##
 ## The optional return value @var{h} is a graphics handle to the created plot.
 ##
+## Example:
+##
 ## @example
 ## @group
 ## f = @@(x,y) sqrt (abs (x .* y)) ./ (1 + x.^2 + y.^2);
@@ -44,19 +49,19 @@
 ## @end group
 ## @end example
 ##
-## @seealso{ezplot, ezcontour, ezsurfc, ezmeshc}
+## @seealso{contourf, ezcontour, ezplot, ezmeshc, ezsurfc}
 ## @end deftypefn
 
-function retval = ezcontourf (varargin)
+function h = ezcontourf (varargin)
 
-  [h, needusage] = __ezplot__ ("contourf", varargin{:});
+  [htmp, needusage] = __ezplot__ ("contourf", varargin{:});
 
   if (needusage)
     print_usage ();
   endif
 
   if (nargout > 0)
-    retval = h;
+    h = htmp;
   endif
 
 endfunction

--- a/scripts/plot/ezmesh.m	Tue Jun 04 19:41:21 2013 +0200
+++ b/scripts/plot/ezmesh.m	Wed Jul 17 10:09:44 2013 -0700
@@ -22,31 +22,36 @@
 ## @deftypefnx {Function File} {} ezmesh (@dots{}, @var{dom})
 ## @deftypefnx {Function File} {} ezmesh (@dots{}, @var{n})
 ## @deftypefnx {Function File} {} ezmesh (@dots{}, "circ")
-## @deftypefnx {Function File} {} ezmesh (@var{h}, @dots{})
+## @deftypefnx {Function File} {} ezmesh (@var{hax}, @dots{})
 ## @deftypefnx {Function File} {@var{h} =} ezmesh (@dots{})
 ##
-## Plot the mesh defined by a function.  @var{f} is a string, inline
-## function or function handle with two arguments defining the function.  By
-## default the plot is over the domain @code{-2*pi < @var{x} < 2*pi} and
-## @code{-2*pi < @var{y} < 2*pi} with 60 points in each dimension.
+## Plot the mesh defined by a function.
 ##
-## If @var{dom} is a two element vector, it represents the minimum and maximum
-## value of both @var{x} and @var{y}.  If @var{dom} is a four element vector,
-## then the minimum and maximum value of @var{x} and @var{y} are specify
-## separately.
-##
-## @var{n} is a scalar defining the number of points to use in each dimension.
+## @var{f} is a string, inline function, or function handle with two arguments
+## defining the function.  By default the plot is over the meshed domain
+## @code{-2*pi <= @var{x} | @var{y} <= 2*pi} with 60 points in each dimension.
 ##
 ## If three functions are passed, then plot the parametrically defined
 ## function @code{[@var{fx} (@var{s}, @var{t}), @var{fy} (@var{s}, @var{t}),
 ## @var{fz} (@var{s}, @var{t})]}.
 ##
+## If @var{dom} is a two element vector, it represents the minimum and maximum
+## values of both @var{x} and @var{y}.  If @var{dom} is a four element vector,
+## then the minimum and maximum values are @code{[xmin xmax ymin ymax]}.
+##
+## @var{n} is a scalar defining the number of points to use in each dimension.
+##
 ## If the argument "circ" is given, then the function is plotted over a disk
 ## centered on the middle of the domain @var{dom}.
 ##
+## If the first argument is an axis handle, @var{hax}, then plot into this
+## axis rather than the current axis handle returned by @code{gca}.
+##
 ## The optional return value @var{h} is a graphics handle to the created 
 ## surface object.
 ##
+## Example 1: 2-argument function
+##
 ## @example
 ## @group
 ## f = @@(x,y) sqrt (abs (x .* y)) ./ (1 + x.^2 + y.^2);
@@ -54,7 +59,7 @@
 ## @end group
 ## @end example
 ##
-## An example of a parametrically defined function is
+## Example 2: parametrically defined function
 ##
 ## @example
 ## @group
@@ -65,19 +70,19 @@
 ## @end group
 ## @end example
 ##
-## @seealso{ezplot, ezmeshc, ezsurf, ezsurfc}
+## @seealso{mesh, ezmeshc, ezplot, ezsurf, ezsurfc, hidden}
 ## @end deftypefn
 
-function retval = ezmesh (varargin)
+function h = ezmesh (varargin)
 
-  [h, needusage] = __ezplot__ ("mesh", varargin{:});
+  [htmp, needusage] = __ezplot__ ("mesh", varargin{:});
 
   if (needusage)
     print_usage ();
   endif
 
   if (nargout > 0)
-    retval = h;
+    h = htmp;
   endif
 
 endfunction

--- a/scripts/plot/ezmeshc.m	Tue Jun 04 19:41:21 2013 +0200
+++ b/scripts/plot/ezmeshc.m	Wed Jul 17 10:09:44 2013 -0700
@@ -22,25 +22,25 @@
 ## @deftypefnx {Function File} {} ezmeshc (@dots{}, @var{dom})
 ## @deftypefnx {Function File} {} ezmeshc (@dots{}, @var{n})
 ## @deftypefnx {Function File} {} ezmeshc (@dots{}, "circ")
-## @deftypefnx {Function File} {} ezmeshc (@var{h}, @dots{})
+## @deftypefnx {Function File} {} ezmeshc (@var{hax}, @dots{})
 ## @deftypefnx {Function File} {@var{h} =} ezmeshc (@dots{})
 ##
-## Plot the mesh and contour lines defined by a function.  @var{f} is a string,
-## inline function or function handle with two arguments defining the function.
-## By default the plot is over the domain @code{-2*pi < @var{x} < 2*pi} and
-## @code{-2*pi < @var{y} < 2*pi} with 60 points in each dimension.
+## Plot the mesh and contour lines defined by a function.
 ##
-## If @var{dom} is a two element vector, it represents the minimum and maximum
-## value of both @var{x} and @var{y}.  If @var{dom} is a four element vector,
-## then the minimum and maximum value of @var{x} and @var{y} are specify
-## separately.
-##
-## @var{n} is a scalar defining the number of points to use in each dimension.
+## @var{f} is a string, inline function, or function handle with two arguments
+## defining the function.  By default the plot is over the meshed domain
+## @code{-2*pi <= @var{x} | @var{y} <= 2*pi} with 60 points in each dimension.
 ##
 ## If three functions are passed, then plot the parametrically defined
 ## function @code{[@var{fx} (@var{s}, @var{t}), @var{fy} (@var{s}, @var{t}),
 ## @var{fz} (@var{s}, @var{t})]}.
 ##
+## If @var{dom} is a two element vector, it represents the minimum and maximum
+## values of both @var{x} and @var{y}.  If @var{dom} is a four element vector,
+## then the minimum and maximum values are @code{[xmin xmax ymin ymax]}.
+##
+## @var{n} is a scalar defining the number of points to use in each dimension.
+##
 ## If the argument "circ" is given, then the function is plotted over a disk
 ## centered on the middle of the domain @var{dom}.
 ##
@@ -48,6 +48,8 @@
 ## handle for the created mesh plot and a second handle for the created contour
 ## plot.
 ##
+## Example: 2-argument function
+##
 ## @example
 ## @group
 ## f = @@(x,y) sqrt (abs (x .* y)) ./ (1 + x.^2 + y.^2);
@@ -55,19 +57,19 @@
 ## @end group
 ## @end example
 ##
-## @seealso{ezplot, ezsurfc, ezsurf, ezmesh}
+## @seealso{meshc, ezmesh, ezplot, ezsurf, ezsurfc, hidden}
 ## @end deftypefn
 
-function retval = ezmeshc (varargin)
+function h = ezmeshc (varargin)
 
-  [h, needusage] = __ezplot__ ("meshc", varargin{:});
+  [htmp, needusage] = __ezplot__ ("meshc", varargin{:});
 
   if (needusage)
     print_usage ();
   endif
 
   if (nargout > 0)
-    retval = h;
+    h = htmp;
   endif
 
 endfunction

--- a/scripts/plot/ezplot.m	Tue Jun 04 19:41:21 2013 +0200
+++ b/scripts/plot/ezplot.m	Wed Jul 17 10:09:44 2013 -0700
@@ -18,28 +18,31 @@
 
 ## -*- texinfo -*-
 ## @deftypefn  {Function File} {} ezplot (@var{f})
+## @deftypefnx {Function File} {} ezplot (@var{f2v})
 ## @deftypefnx {Function File} {} ezplot (@var{fx}, @var{fy})
 ## @deftypefnx {Function File} {} ezplot (@dots{}, @var{dom})
 ## @deftypefnx {Function File} {} ezplot (@dots{}, @var{n})
-## @deftypefnx {Function File} {} ezplot (@var{h}, @dots{})
+## @deftypefnx {Function File} {} ezplot (@var{hax}, @dots{})
 ## @deftypefnx {Function File} {@var{h} =} ezplot (@dots{})
 ##
-## Plot the curve defined by @var{f} in two dimensions.  The function
-## @var{f} may be a string, inline function or function handle and can
-## have either one or two variables.  If @var{f} has one variable, then
+## Plot the 2-D curve defined by the function @var{f}.
+##
+## The function @var{f} may be a string, inline function, or function handle
+## and can have either one or two variables.  If @var{f} has one variable, then
 ## the function is plotted over the domain @code{-2*pi < @var{x} < 2*pi}
 ## with 500 points.
 ##
-## If @var{f} has two variables then @code{@var{f}(@var{x},@var{y}) = 0}
-## is calculated over the meshed domain @code{-2*pi < @var{x} | @var{y}
-## < 2*pi} with 60 by 60 in the mesh.  For example:
+## If @var{f2v} is a function of two variables then the implicit function
+## @code{@var{f}(@var{x},@var{y}) = 0} is calculated over the meshed domain
+## @code{-2*pi <= @var{x} | @var{y} <= 2*pi} with 60 points in each dimension.
+##
+## For example:
 ##
 ## @example
 ## ezplot (@@(@var{x}, @var{y}) @var{x}.^2 - @var{y}.^2 - 1)
 ## @end example
 ##
-## If two functions are passed as strings, inline functions or function
-## handles, then the parametric function
+## If two functions are passed as inputs then the parametric function
 ##
 ## @example
 ## @group
@@ -49,47 +52,57 @@
 ## @end example
 ##
 ## @noindent
-## is plotted over the domain @code{-2*pi < @var{t} < 2*pi} with 500
-## points.
+## is plotted over the domain @code{-2*pi <= @var{t} <= 2*pi} with 500 points.
 ##
 ## If @var{dom} is a two element vector, it represents the minimum and maximum
-## value of @var{x}, @var{y} and @var{t}.  If it is a four element
-## vector, then the minimum and maximum values of @var{x} and @var{t}
-## are determined by the first two elements and the minimum and maximum
-## of @var{y} by the second pair of elements.
+## values of both @var{x} and @var{y}, or @var{t} for a parametric plot.  If
+## @var{dom} is a four element vector, then the minimum and maximum values are
+## @code{[xmin xmax ymin ymax]}.
 ##
 ## @var{n} is a scalar defining the number of points to use in plotting
 ## the function.
 ##
-## The optional return value @var{h} is a graphics handle to the created plot.
+## If the first argument is an axis handle, @var{hax}, then plot into this
+## axis rather than the current axis handle returned by @code{gca}.
 ##
-## @seealso{plot, ezplot3}
+## The optional return value @var{h} is a list of graphics handles in the
+## created plot.
+##
+## @seealso{plot, ezplot3, ezpolar, ezcontour, ezcontourf, ezmesh, ezmeshc, ezsurf, ezsurfc}
 ## @end deftypefn
 
-function retval = ezplot (varargin)
+function h = ezplot (varargin)
 
-  [h, needusage] = __ezplot__ ("plot", varargin{:});
+  [htmp, needusage] = __ezplot__ ("plot", varargin{:});
 
   if (needusage)
     print_usage ();
   endif
 
   if (nargout > 0)
-    retval = h;
+    h = htmp;
   endif
 
 endfunction
 
 
 %!demo
+%! ## sinc function using function handle
+%! f = @(x) sin (pi*x) ./ (pi*x);
+%! ezplot (f);
+
+%!demo
+%! ## example of a function string and explicit limits
+%! clf;
+%! ezplot ('1/x', [-2 2]);
+
+%!demo
+%! ## parameterized function example over -2*pi <= t <= +2*pi
 %! clf;
 %! ezplot (@cos, @sin);
 
 %!demo
+%! ## implicit function of 2 variables
 %! clf;
-%! ezplot ('1/x');
+%! ezplot (inline ('x^2 - y^2 - 1'));
 
-%!demo
-%! clf;
-%! ezplot (inline ('x^2 - y^2 = 1'));
-

--- a/scripts/plot/ezplot3.m	Tue Jun 04 19:41:21 2013 +0200
+++ b/scripts/plot/ezplot3.m	Wed Jul 17 10:09:44 2013 -0700
@@ -20,17 +20,24 @@
 ## @deftypefn  {Function File} {} ezplot3 (@var{fx}, @var{fy}, @var{fz})
 ## @deftypefnx {Function File} {} ezplot3 (@dots{}, @var{dom})
 ## @deftypefnx {Function File} {} ezplot3 (@dots{}, @var{n})
-## @deftypefnx {Function File} {} ezplot3 (@var{h}, @dots{})
+## @deftypefnx {Function File} {} ezplot3 (@var{hax}, @dots{})
 ## @deftypefnx {Function File} {@var{h} =} ezplot3 (@dots{})
 ##
 ## Plot a parametrically defined curve in three dimensions.
-## @var{fx}, @var{fy}, and @var{fz} are strings, inline functions
-## or function handles with one arguments defining the function.  By
-## default the plot is over the domain @code{-2*pi < @var{x} < 2*pi}
-## with 60 points.
+##
+## @var{fx}, @var{fy}, and @var{fz} are strings, inline functions,
+## or function handles with one argument defining the function.  By
+## default the plot is over the domain @code{0 <= @var{t} <= 2*pi}
+## with 500 points.
 ##
 ## If @var{dom} is a two element vector, it represents the minimum and maximum
-## value of @var{t}.  @var{n} is a scalar defining the number of points to use.
+## values of @var{t}.
+##
+## @var{n} is a scalar defining the number of points to use in plotting the
+## function.
+##
+## If the first argument is an axis handle, @var{hax}, then plot into this
+## axis rather than the current axis handle returned by @code{gca}.
 ##
 ## The optional return value @var{h} is a graphics handle to the created plot.
 ##
@@ -43,19 +50,19 @@
 ## @end group
 ## @end example
 ##
-## @seealso{plot3, ezplot, ezsurf, ezmesh}
+## @seealso{plot3, ezplot, ezmesh, ezsurf}
 ## @end deftypefn
 
-function retval = ezplot3 (varargin)
+function h = ezplot3 (varargin)
 
-  [h, needusage] = __ezplot__ ("plot3", varargin{:});
+  [htmp, needusage] = __ezplot__ ("plot3", varargin{:});
 
   if (needusage)
     print_usage ();
   endif
 
   if (nargout > 0)
-    retval = h;
+    h = htmp;
   endif
 
 endfunction
@@ -68,3 +75,10 @@
 %! fz = @(t) t;
 %! ezplot3 (fx, fy, fz, [0, 10*pi], 100);
 
+%!demo
+%! clf;
+%! fx = @(t) cos (t);
+%! fy = @(t) sin (t);
+%! fz = @(t) t;
+%! ezplot3 (fx, fy, fz, [0, 10*pi], 100, 'animate');
+

--- a/scripts/plot/ezpolar.m	Tue Jun 04 19:41:21 2013 +0200
+++ b/scripts/plot/ezpolar.m	Wed Jul 17 10:09:44 2013 -0700
@@ -23,17 +23,22 @@
 ## @deftypefnx {Function File} {} ezpolar (@var{hax}, @dots{})
 ## @deftypefnx {Function File} {@var{h} =} ezpolar (@dots{})
 ##
-## Plot a function in polar coordinates.  The function @var{f} is
-## a string, inline function, or function handle with a single argument.
-## The expected form of the function is
+## Plot a 2-D function in polar coordinates.
+## 
+## The function @var{f} is a string, inline function, or function handle with
+## a single argument.  The expected form of the function is
 ## @code{@var{rho} = @var{f}(@var{theta})}.
-## By default the plot is over the domain @code{0 < @var{theta} < 2*pi} with 60
-## points.
+## By default the plot is over the domain @code{0 <= @var{theta} <= 2*pi}
+## with 500 points.
 ##
 ## If @var{dom} is a two element vector, it represents the minimum and maximum
-## values of @var{theta}.  @var{n} is a scalar defining the number of points to
-## use.  If the optional input @var{hax} is given then the plot is placed into
-## the specified axes rather than the current axes.
+## values of @var{theta}.
+##
+## @var{n} is a scalar defining the number of points to use in plotting
+## the function.
+##
+## If the first argument is an axis handle, @var{hax}, then plot into this
+## axis rather than the current axis handle returned by @code{gca}.
 ##
 ## The optional return value @var{h} is a graphics handle to the created plot.
 ##
@@ -43,19 +48,19 @@
 ## ezpolar (@@(t) 1 + sin (t));
 ## @end example
 ##
-## @seealso{polar, ezplot, ezsurf, ezmesh}
+## @seealso{polar, ezplot}
 ## @end deftypefn
 
-function retval = ezpolar (varargin)
+function h = ezpolar (varargin)
 
-  [h, needusage] = __ezplot__ ("polar", varargin{:});
+  [htmp, needusage] = __ezplot__ ("polar", varargin{:});
 
   if (needusage)
     print_usage ();
   endif
 
   if (nargout > 0)
-    retval = h;
+    h = htmp;
   endif
 
 endfunction

--- a/scripts/plot/ezsurf.m	Tue Jun 04 19:41:21 2013 +0200
+++ b/scripts/plot/ezsurf.m	Wed Jul 17 10:09:44 2013 -0700
@@ -22,31 +22,36 @@
 ## @deftypefnx {Function File} {} ezsurf (@dots{}, @var{dom})
 ## @deftypefnx {Function File} {} ezsurf (@dots{}, @var{n})
 ## @deftypefnx {Function File} {} ezsurf (@dots{}, "circ")
-## @deftypefnx {Function File} {} ezsurf (@var{h}, @dots{})
+## @deftypefnx {Function File} {} ezsurf (@var{hax}, @dots{})
 ## @deftypefnx {Function File} {@var{h} =} ezsurf (@dots{})
 ##
-## Plot the surface defined by a function.  @var{f} is a string, inline
-## function or function handle with two arguments defining the function.  By
-## default the plot is over the domain @code{-2*pi < @var{x} < 2*pi} and
-## @code{-2*pi < @var{y} < 2*pi} with 60 points in each dimension.
+## Plot the surface defined by a function.
 ##
-## If @var{dom} is a two element vector, it represents the minimum and maximum
-## value of both @var{x} and @var{y}.  If @var{dom} is a four element vector,
-## then the minimum and maximum value of @var{x} and @var{y} are specify
-## separately.
-##
-## @var{n} is a scalar defining the number of points to use in each dimension.
+## @var{f} is a string, inline function, or function handle with two arguments
+## defining the function.  By default the plot is over the meshed domain
+## @code{-2*pi <= @var{x} | @var{y} <= 2*pi} with 60 points in each dimension.
 ##
 ## If three functions are passed, then plot the parametrically defined
 ## function @code{[@var{fx} (@var{s}, @var{t}), @var{fy} (@var{s}, @var{t}),
 ## @var{fz} (@var{s}, @var{t})]}.
 ##
+## If @var{dom} is a two element vector, it represents the minimum and maximum
+## values of both @var{x} and @var{y}.  If @var{dom} is a four element vector,
+## then the minimum and maximum values are @code{[xmin xmax ymin ymax]}.
+##
+## @var{n} is a scalar defining the number of points to use in each dimension.
+##
 ## If the argument "circ" is given, then the function is plotted over a disk
 ## centered on the middle of the domain @var{dom}.
 ##
+## If the first argument is an axis handle, @var{hax}, then plot into this
+## axis rather than the current axis handle returned by @code{gca}.
+##
 ## The optional return value @var{h} is a graphics handle to the created
 ## surface object.
 ##
+## Example 1: 2-argument function
+##
 ## @example
 ## @group
 ## f = @@(x,y) sqrt (abs (x .* y)) ./ (1 + x.^2 + y.^2);
@@ -54,7 +59,7 @@
 ## @end group
 ## @end example
 ##
-## An example of a parametrically defined function is
+## Example 2: parametrically defined function
 ##
 ## @example
 ## @group
@@ -65,19 +70,19 @@
 ## @end group
 ## @end example
 ##
-## @seealso{ezplot, ezmesh, ezsurfc, ezmeshc}
+## @seealso{surf, ezsurfc, ezplot, ezmesh, ezmeshc, shading}
 ## @end deftypefn
 
-function retval = ezsurf (varargin)
+function h = ezsurf (varargin)
 
-  [h, needusage] = __ezplot__ ("surf", varargin{:});
+  [htmp, needusage] = __ezplot__ ("surf", varargin{:});
 
   if (needusage)
     print_usage ();
   endif
 
   if (nargout > 0)
-    retval = h;
+    h = htmp;
   endif
 
 endfunction
@@ -96,4 +101,16 @@
 %! fy = @(s,t) sin (s) .* cos (t);
 %! fz = @(s,t) sin (t);
 %! ezsurf (fx, fy, fz, [-pi,pi,-pi/2,pi/2], 20);
+%! axis equal;
 
+%!demo
+%! clf;
+%! colormap ('default');
+%! f = @(x,y) x.^2 + y.^2;
+%! subplot (1,2,1);
+%!  ezsurf (f, [-2,2]);
+%!  title ({'x^2 + y^2'; 'plotted over rectangular grid (default)'});
+%! subplot (1,2,2);
+%!  ezsurf (f, [-2,2], 'circ');
+%!  title ({'x^2 + y^2'; 'plotted over circular disk with "circ"'});
+

--- a/scripts/plot/ezsurfc.m	Tue Jun 04 19:41:21 2013 +0200
+++ b/scripts/plot/ezsurfc.m	Wed Jul 17 10:09:44 2013 -0700
@@ -22,31 +22,36 @@
 ## @deftypefnx {Function File} {} ezsurfc (@dots{}, @var{dom})
 ## @deftypefnx {Function File} {} ezsurfc (@dots{}, @var{n})
 ## @deftypefnx {Function File} {} ezsurfc (@dots{}, "circ")
-## @deftypefnx {Function File} {} ezsurfc (@var{h}, @dots{})
+## @deftypefnx {Function File} {} ezsurfc (@var{hax}, @dots{})
 ## @deftypefnx {Function File} {@var{h} =} ezsurfc (@dots{})
 ##
-## Plot the surface and contour lines defined by a function.  @var{f} is a
-## string, inline function or function handle with two arguments defining the
-## function.  By default the plot is over the domain @code{-2*pi < @var{x} <
-## 2*pi} and @code{-2*pi < @var{y} < 2*pi} with 60 points in each dimension.
+## Plot the surface and contour lines defined by a function.
 ##
-## If @var{dom} is a two element vector, it represents the minimum and maximum
-## value of both @var{x} and @var{y}.  If @var{dom} is a four element vector,
-## then the minimum and maximum value of @var{x} and @var{y} are specify
-## separately.
-##
-## @var{n} is a scalar defining the number of points to use in each dimension.
+## @var{f} is a string, inline function, or function handle with two arguments
+## defining the function.  By default the plot is over the meshed domain
+## @code{-2*pi <= @var{x} | @var{y} <= 2*pi} with 60 points in each dimension.
 ##
 ## If three functions are passed, then plot the parametrically defined
 ## function @code{[@var{fx} (@var{s}, @var{t}), @var{fy} (@var{s}, @var{t}),
 ## @var{fz} (@var{s}, @var{t})]}.
 ##
+## If @var{dom} is a two element vector, it represents the minimum and maximum
+## values of both @var{x} and @var{y}.  If @var{dom} is a four element vector,
+## then the minimum and maximum values are @code{[xmin xmax ymin ymax]}.
+##
+## @var{n} is a scalar defining the number of points to use in each dimension.
+##
 ## If the argument "circ" is given, then the function is plotted over a disk
 ## centered on the middle of the domain @var{dom}.
 ##
+## If the first argument is an axis handle, @var{hax}, then plot into this
+## axis rather than the current axis handle returned by @code{gca}.
+##
 ## The optional return value @var{h} is a 2-element vector with a graphics
-## for the created surface plot and a second handle for the created contour
-## plot.
+## handle for the created surface plot and a second handle for the created
+## contour plot.
+##
+## Example:
 ##
 ## @example
 ## @group
@@ -55,19 +60,19 @@
 ## @end group
 ## @end example
 ##
-## @seealso{ezplot, ezmeshc, ezsurf, ezmesh}
+## @seealso{surfc, ezsurf, ezplot, ezmesh, ezmeshc, shading}
 ## @end deftypefn
 
-function retval = ezsurfc (varargin)
+function h = ezsurfc (varargin)
 
-  [h, needusage] = __ezplot__ ("surfc", varargin{:});
+  [htmp, needusage] = __ezplot__ ("surfc", varargin{:});
 
   if (needusage)
     print_usage ();
   endif
 
   if (nargout > 0)
-    retval = h;
+    h = htmp;
   endif
 
 endfunction

--- a/scripts/plot/private/__ezplot__.m	Tue Jun 04 19:41:21 2013 +0200
+++ b/scripts/plot/private/__ezplot__.m	Wed Jul 17 10:09:44 2013 -0700
@@ -17,66 +17,75 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{h}, @var{needusage}] =} __ezplot__ (@var{pfunc}, @var{varargin})
+## @deftypefn {Function File} {[@var{h}, @var{needusage}] =} __ezplot__ (@var{pltfunc}, @var{varargin})
 ## Undocumented internal function.
 ## @end deftypefn
 
-function [h, needusage] = __ezplot__ (pfunc, varargin)
+## Overview: This function is the back-end for the 9 ez* plot functions.
+##           As such, most of the function is actually dedicated to sorting
+##           out the inputs and verifying that the particular ez* function
+##           called was called correctly.  The actual plotting occurs near
+##           the end in an unwind_protect block. 
 
-  func = cstrcat ("ez", pfunc);
-  if (strncmp (pfunc, "contour", 7))
-    iscontour = true;
-  else
-    iscontour = false;
-  endif
-  if (strcmp (pfunc, "plot"))
-    isplot = true;
-    isplot3 = false;
-    ispolar = false;
-    nargs = 1;
-  elseif (strcmp (pfunc, "plot3"))
-    isplot = false;
-    isplot3 = true;
-    ispolar = false;
-    nargs = 1;
-  elseif (strcmp (pfunc, "polar"))
-    isplot = false;
-    isplot3 = false;
-    ispolar = true;
-    nargs = 1;
-  else
-    isplot = false;
-    isplot3 = false;
-    ispolar = false;
-    nargs = 2;
-  endif
+function [h, needusage] = __ezplot__ (pltfunc, varargin)
+
+  ezfunc = ["ez" pltfunc];
 
-  [ax, varargin, nargin] = __plt_get_axis_arg__ (func, varargin{:});
+  [ax, varargin, nargin] = __plt_get_axis_arg__ (ezfunc, varargin{:});
 
+  ## Define outputs early in case of shorting out of function with return;
+  h = [];
   needusage = false;
   if (nargin < 1)
     needusage = true;
     return;
   endif
 
+  iscontour = strncmp (pltfunc, "contour", 7);
+
+  ## Defaults for ezplot
+  isplot  = true;
+  isplot3 = false;
+  ispolar = false;
+  nargs = 1;
+  switch (pltfunc)
+    case "plot"
+      ## defaults already set
+
+    case "plot3"
+      isplot  = false;
+      isplot3 = true;
+  
+    case "polar"
+      isplot  = false;
+      ispolar = true;
+    
+    otherwise
+      ## contour, mesh, surf plots
+      isplot  = false;
+      nargs = 2;
+
+  endswitch
+
   parametric = false;
-  fun = varargin {1};
+  fun = varargin{1};
   if (ischar (fun))
     if (exist (fun, "file") || exist (fun, "builtin"))
-      fun = vectorize (inline (cstrcat (fun, "(t)")));
+      fun = inline ([fun "(t)"]);
     else
       fun = vectorize (inline (fun));
     endif
-    if (isplot && length (argnames (fun)) == 2)
+    argids = argnames (fun);
+    if (isplot && length (argids) == 2)
       nargs = 2;
-    elseif (length (argnames (fun)) != nargs)
-      error ("%s: excepting a function of %d arguments", func, nargs);
+    elseif (numel (argids) != nargs)
+      error ("%s: expecting a function of %d arguments", ezfunc, nargs);
     endif
     fstr = formula (fun);
     if (isplot)
-      xarg = (argnames (fun)){1};
+      xarg = argids{1};
       if (nargs == 2)
-        yarg = (argnames (fun)){2};
+        yarg = argids{2};
       else
         yarg = "";
       endif
@@ -87,21 +96,22 @@
       xarg = "";
       yarg = "";
     else
-      xarg = (argnames (fun)){1};
-      yarg = (argnames (fun)){2};
+      xarg = argids{1};
+      yarg = argids{2};
     endif
   elseif (strcmp (typeinfo (fun), "inline function"))
-    if (isplot && length (argnames (fun)) == 2)
+    argids = argnames (fun);
+    if (isplot && length (argids) == 2)
       nargs = 2;
-    elseif (length (argnames (fun)) != nargs)
-      error ("%s: excepting a function of %d arguments", func, nargs);
+    elseif (numel (argids) != nargs)
+      error ("%s: expecting a function of %d arguments", ezfunc, nargs);
     endif
     fun = vectorize (fun);
     fstr = formula (fun);
     if (isplot)
-      xarg = (argnames (fun)){1};
+      xarg = argids{1};
       if (nargs == 2)
-        yarg = (argnames (fun)){2};
+        yarg = argids{2};
       else
         yarg = "";
       endif
@@ -112,27 +122,27 @@
       xarg = "";
       yarg = "";
     else
-      xarg = (argnames (fun))(1);
-      yarg = (argnames (fun))(2);
+      xarg = argids{1};
+      yarg = argids{2};
     endif
   elseif (isa (fun, "function_handle"))
     fstr = func2str (fun);
-    if (! isempty (strfind (fstr, ')')))
-      args = regexp (substr (fstr, 3, strfind (fstr, ')')(1) - 3),
-                     '(\w+)', 'tokens');
-    fstr = substr (fstr, strfind (fstr, ')')(1) + 1);
+    idx = index (fstr, ')');
+    if (idx != 0)
+      args = regexp (fstr(3:(idx-1)), '\w+', 'match');
+      fstr = fstr(idx+2:end);  # remove '@(x) ' from string name
     else
-      args = {{"x"}};
+      args = {"x"};
     endif
     if (isplot && length (args) == 2)
       nargs = 2;
-    elseif (length (args) != nargs)
-      error ("%s: excepting a function of %d arguments", func, nargs);
+    elseif (numel (args) != nargs)
+      error ("%s: expecting a function of %d arguments", ezfunc, nargs);
     endif
     if (isplot)
-      xarg = args{1}{1};
+      xarg = args{1};
       if (nargs == 2)
-        yarg = args{2}{1};
+        yarg = args{2};
       else
         yarg = "";
       endif
@@ -143,95 +153,101 @@
       xarg = "";
       yarg = "";
     else
-      xarg = args{1}{1};
-      yarg = args{2}{1};
+      xarg = args{1};
+      yarg = args{2};
     endif
   else
-    error ("%s: expecting string, inline function or function handle", func);
+    error ("%s: expecting string, inline function, or function handle", ezfunc);
   endif
 
   if (nargin > 2 || (nargin == 2 && isplot))
     funx = fun;
     fstrx = fstr;
-    funy = varargin {2};
+    funy = varargin{2};
     if (ischar (funy) && ! strcmp (funy, "circ") && ! strcmp (funy, "animate"))
       parametric = true;
       if (exist (funy, "file") || exist (funy, "builtin"))
-        funy = vectorize (inline (cstrcat (funy, "(t)")));
+        funy = inline ([funy "(t)"]);
       else
         funy = vectorize (inline (funy));
       endif
-      if (length (argnames (funy)) != nargs)
-        error ("%s: excepting a function of %d arguments", func, nargs);
+      if (numel (argnames (funy)) != nargs)
+        error ("%s: expecting a function of %d arguments", ezfunc, nargs);
       endif
       fstry = formula (funy);
     elseif (strcmp (typeinfo (funy), "inline function"))
       parametric = true;
-      if (length (argnames (funy)) != nargs)
-        error ("%s: excepting a function of %d arguments", func, nargs);
+      if (numel (argnames (funy)) != nargs)
+        error ("%s: expecting a function of %d arguments", ezfunc, nargs);
       endif
       funy = vectorize (funy);
       fstry = formula (funy);
     elseif (isa (funy, "function_handle"))
       parametric = true;
       fstry = func2str (funy);
-      if (! isempty (strfind (fstry, ')')))
-        args = regexp (substr (fstry, 3, strfind (fstry, ')')(1) - 3),
-                       '(\w+)', 'tokens');
-        fstry = substr (fstry, strfind (fstry, ')')(1) + 1);
+      idx = index (fstry, ')');
+      if (idx != 0)
+        args = regexp (fstry(3:(idx-1)), '\w+', 'match');
+        fstry = fstry(idx+2:end);  # remove '@(x) ' from string name
       else
-        args = {{"y"}};
+        args = {"y"};
       endif
-      if (length (args) != nargs)
-        error ("%s: excepting a function of %d arguments", func, nargs);
+      if (numel (args) != nargs)
+        error ("%s: expecting a function of %d arguments", ezfunc, nargs);
       endif
     endif
 
-    if (parametric && isplot)
-      xarg = "x";
-      yarg = "y";
+    if (! parametric && isplot3)
+      needusage = true;  # Can't call non-parametric ezplot3
+      return;
+    elseif (parametric && isplot)
       if (nargs == 2)
-        error ("%s: can not define a parametric function in this manner");
+        error ("%s: can not define a parametric function in this manner", ezfunc);
+      else
+        xarg = "x";
+        yarg = "y";
       endif
-    endif
-
-    if (!isplot && parametric)
-      funz = varargin {3};
+    elseif (parametric)
+      funz = varargin{3};
       if (ischar (funz) && ! strcmp (funz, "circ")
           && ! strcmp (funz, "animate"))
         if (exist (funz, "file") || exist (funz, "builtin"))
-          funz = vectorize (inline (cstrcat (funz, "(t)")));
+          funz = inline ([funz "(t)"]);
         else
           funz = vectorize (inline (funz));
         endif
-        if (length (argnames (funz)) != nargs)
-          error ("%s: excepting a function of %d arguments", func, nargs);
+        if (numel (argnames (funz)) > nargs)
+          error ("%s: expecting a function of %d arguments", ezfunc, nargs);
         endif
         fstrz = formula (funz);
       elseif (strcmp (typeinfo (funz), "inline function"))
-        if (length (argnames (funz)) != nargs)
-          error ("%s: excepting a function of %d arguments", func, nargs);
+        if (numel (argnames (funz)) != nargs)
+          error ("%s: expecting a function of %d arguments", ezfunc, nargs);
         endif
         funz = vectorize (funz);
         fstrz = formula (funz);
       elseif (isa (funz, "function_handle"))
         fstrz = func2str (funz);
-        args = regexp (substr (fstrz, 3, strfind (fstrz, ')')(1) - 3),
-                       '(\w+)', 'tokens');
-        if (length (args) != nargs)
-          error ("%s: excepting a function of %d arguments", func, nargs);
+        idx = index (fstrz, ')');
+        if (idx != 0)
+          args = regexp (fstrz(3:(idx-1)), '\w+', 'match');
+          fstrz = fstrz(idx+2:end);  # remove '@(x) ' from string name
+        else
+          args = {"z"};
         endif
-        fstrz = substr (fstrz, strfind (fstrz, ')')(1) + 1);
+        if (numel (args) != nargs)
+          error ("%s: expecting a function of %d arguments", ezfunc, nargs);
+        endif
       else
-        error ("%s: parametric plots expect 3 functions", func);
+        error ("%s: parametric plots expect 3 functions", ezfunc);
       endif
     endif
   endif
 
-  if (isplot && nargs != 2)
-    n = 500;
+  if ((isplot && nargs != 2) || isplot3 || ispolar)
+    n = 500;   # default for point-style functions like plot
   else
-    n = 60;
+    n = 60;    # default for meshgrid style functions like contour, surf
   endif
   domain = [];
   circ = false;
@@ -254,181 +270,201 @@
     elseif (isscalar (arg))
       n = arg;
     elseif (numel (arg) == 2)
-      domain = [arg(:).' arg(:).'];
+      domain = [arg(1) arg(2) arg(1) arg(2)];
     elseif (numel (arg) == 4)
       domain = arg(:).';
     else
-      error ("%s: expecting scalar, 2 or 4 element vector", func);
+      error ("%s: expecting scalar, 2-, or 4-element vector", ezfunc);
     endif
   endwhile
 
+  if (circ && (iscontour || isplot3 || isplot))
+    needusage = true;
+    return;
+  elseif (circ && parametric)
+    error ("%s: can not have both circular domain and parametric function",
+           ezfunc);
+  endif
+
+  if (animate && ! isplot3)
+    error ("%s: animate option only valid for ezplot3", ezfunc);
+  endif
+
+  if (parametric)
+    ## Make the label strings pretty by removing extra spaces between base
+    ## and exponent, the '.' in vectorized code, and the '*' for multiply.
+    fstrx = regexprep (regexprep (regexprep (fstrx,
+           '\s*\.?(?:\^|\*\*)\s*','^'), '\.([/+-])', '$1'), '\s*\.?\*\s*', ' ');
+    fstry = regexprep (regexprep (regexprep (fstry,
+           '\s*\.?(?:\^|\*\*)\s*','^'), '\.([/+-])', '$1'), '\s*\.?\*\s*', ' ');
+    if (isplot)
+      fstr = ["x = " fstrx ", y = " fstry];
+    else
+      fstrz = regexprep (regexprep (regexprep (fstrz,
+           '\s*\.?(?:\^|\*\*)\s*','^'), '\.([/+-])', '$1'), '\s*\.?\*\s*', ' ');
+      fstr = ["x = " fstrx ",y = " fstry ", z = " fstrz];
+    endif
+  else
+    fstr = regexprep (regexprep (regexprep (fstr,
+           '\s*\.?(?:\^|\*\*)\s*','^'), '\.([/+-])', '$1'), '\s*\.?\*\s*', ' ');
+    if (isplot && nargs == 2)
+      fstr = [fstr " = 0"];  # make title string of implicit function
+    endif
+  endif
+
   if (isempty (domain))
+    auto_domain = true;
     if (isplot3 || ispolar)
       domain = [0, 2*pi, 0, 2*pi];
     else
       domain = [-2*pi, 2*pi, -2*pi, 2*pi];
     endif
-  endif
-
-  if (circ)
-    if (iscontour || isplot3 || isplot)
-      needusage = true;
-      return;
-    endif
-    if (parametric)
-      error ("%s: can not have both circular domain and parametric function",
-             func);
-    endif
-    cent = [domain(1) + domain(2), domain(3) + domain(4)] / 2;
-    funx = @(r,t) r .* cos (t) + cent (1);
-    funy = @(r,t) r .* sin (t) + cent (2);
-    domain = [0, sqrt((domain(2) - cent(1))^2 + (domain(4) - cent(2))^2), ...
-              -pi, pi];
-    funz = fun;
-    parametric = true;
-  endif
-
-  if (animate)
-    if (!isplot3)
-      error ("%s: animated graphs only valid with plot3", func);
-    endif
-    error ("%s: animated graphs not implemented", func);
+  else
+    auto_domain = false;
   endif
 
-  if (isplot3 || ispolar || (isplot && nargs == 1))
-    X = linspace (domain (1), domain (2), n);
-  elseif (isplot && numel (domain) == 2)
-    x = linspace (domain (1), domain (2), n);
-    [X, Y] = meshgrid (x, x);
-  else
-    x = linspace (domain (1), domain (2), n);
-    y = linspace (domain (3), domain (4), n);
-    [X, Y] = meshgrid (x, y);
-  endif
+  auto_domain_done = false;
+  do
+    domain_ok = true;
 
-  if (parametric)
-    if (isplot)
-      XX = feval (funx, X);
-      Z = feval (funy, X);
-      X = XX;
-    elseif (isplot3)
-      Z = feval (funz, X);
-      XX = feval (funx, X);
-      YY = feval (funy, X);
-      X = XX;
-      Y = YY;
-    else
-      Z = feval (funz, X, Y);
-      XX = feval (funx, X, Y);
-      YY = feval (funy, X, Y);
-      X = XX;
-      Y = YY;
-
-      ## Eliminate the singularities
-      X = __eliminate_sing__ (X);
-      Y = __eliminate_sing__ (Y);
-      Z = __eliminate_sing__ (Z);
+    if ((isplot && nargs == 1) || isplot3 || ispolar)
+      X = linspace (domain(1), domain(2), n);
+    elseif (isplot && numel (domain) == 2)
+      x = linspace (domain(1), domain(2), n);
+      [X, Y] = meshgrid (x, x);
+    elseif (circ)
+      ## To plot on circular domain develop grid in polar coordinates
+      ## and then switch these to Cartesian coordinates.
+      cent = [domain(1) + domain(2), domain(3) + domain(4)] / 2;
+      rmax = sqrt ((domain(2) - cent(1))^2 + (domain(4) - cent(2))^2);
+      r = linspace (0, rmax, n);
+      t = linspace (0, 2*pi, n);
+      [T, R] = meshgrid (t, r);
+      X = R .* cos (T) + cent(1);
+      Y = R .* sin (T) + cent(2);
+      domain = [-rmax+cent(1), +rmax+cent(1), -rmax+cent(2), +rmax+cent(2)];
+    else  # contour, mesh, surf plots
+      x = linspace (domain(1), domain(2), n);
+      y = linspace (domain(3), domain(4), n);
+      [X, Y] = meshgrid (x, y);
     endif
 
-    fstrx = regexprep (regexprep (regexprep (fstrx,'\s*\.?\^\s*','^'),
-                      '\./', '/'), '\.?\*', '');
-    fstry = regexprep (regexprep (regexprep (fstry,'\s*\.?\^\s*','^'),
-                      '\./', '/'), '\.?\*', '');
-    if (isplot)
-      fstr = cstrcat ("x = ",fstrx,", y = ",fstry);
-    else
-      fstrz = regexprep (regexprep (regexprep (fstrz,'\s*\.?\^\s*','^'),
-                                    '\./', '/'), '\.?\*', '');
-      fstr = cstrcat ("x = ",fstrx,",y = ",fstry,", z = ",fstrz);
-    endif
-  else
-    if (isplot3)
-      needusage = true;
-      return;
-    endif
-
-    fstr = regexprep (regexprep (regexprep (fstr,'\s*\.?\^\s*','^'), '\./', '/'),
-                      '\.?\*', '');
-    if (isplot && nargs == 2)
-      if (strcmp (typeinfo (fun), "inline function")
-          && !isempty (strfind (formula (fun) , "=")))
-        fun = inline (cstrcat (strrep (formula (fun), "=", "- ("), ")"));
+    if (parametric)
+      if (isplot)
+        XX = feval (funx, X);
+        Z = feval (funy, X);
+        X = XX;
+      elseif (isplot3)
+        Z = feval (funz, X);
+        XX = feval (funx, X);
+        YY = feval (funy, X);
+        X = XX;
+        Y = YY;
       else
-        fstr = cstrcat (fstr, " = 0");
-      endif
-
-      Z = feval (fun, X, Y);
-
-      ## Matlab returns line objects for this case and so can't call
-      ## contour directly as it returns patch objects to allow colormaps
-      ## to work with contours. Therefore recreate the lines from the
-      ## output for contourc, and store in cell arrays.
-      [c, lev] = contourc (X, Y, Z, [0, 0]);
+        Z = feval (funz, X, Y);
+        XX = feval (funx, X, Y);
+        YY = feval (funy, X, Y);
+        X = XX;
+        Y = YY;
 
-      i1 = 1;
-      XX = {};
-      YY = {};
-      while (i1 < length (c))
-        clev = c(1,i1);
-        clen = c(2,i1);
-        XX = [XX, {c(1, i1+1:i1+clen)}];
-        YY = [YY, {c(2, i1+1:i1+clen)}];
-        i1 += clen+1;
-      endwhile
-    else
-      if (ispolar)
-        Z = feval (fun, X);
-      elseif (isplot)
-        Z = real (feval (fun, X));
-
-        ## Eliminate the singularities. This seems to be what matlab
-        ## does, but can't be sure.
-        XX = sort (Z (isfinite (Z)));
-        if (length (X) > 4)
-          d = XX(fix (7 * length (XX) / 8)) - XX(fix (length (XX) / 8));
-          yrange = [max(XX(1) - d/8, XX(fix (length (XX) / 8)) - d), ...
-                    min(XX(end) + d/8, XX(fix (7 * length (XX) / 8)) + d)];
-        else
-          yrange = [XX(1), XX(end)];
-        endif
-
-        idx = 2 : length (Z);
-        idx = find (((Z(idx) > yrange(2) / 2) & (Z(idx-1) < yrange(1) / 2)) |
-                 ((Z(idx) < yrange(1) / 2) & (Z(idx-1) > yrange (2) / 2)));
-        if (any (idx))
-          Z(idx) = NaN;
-        endif
-      else
+        ## Eliminate the singularities
+        X = __eliminate_sing__ (X);
+        Y = __eliminate_sing__ (Y);
+        Z = __eliminate_sing__ (Z);
+      endif
+    else  ## non-parametric plots
+      if (isplot && nargs == 2)
         Z = feval (fun, X, Y);
 
-        ## Eliminate the singularities
-        Z = __eliminate_sing__ (Z);
+        ## Matlab returns line objects for this case and so can't call
+        ## contour directly as it returns patch objects to allow colormaps
+        ## to work with contours.  Therefore recreate the lines from the
+        ## output for contourc, and store in cell arrays.
+        [c, ~] = contourc (X, Y, Z, [0, 0]);
+
+        i = 1;
+        XX = YY = {};
+        while (i < length (c))
+          clev = c(1,i);
+          clen = c(2,i);
+          XX = [XX, {c(1, i+1:i+clen)}];
+          YY = [YY, {c(2, i+1:i+clen)}];
+          i += clen+1;
+        endwhile
+      else
+        if (ispolar)
+          Z = feval (fun, X);
+          ## FIXME: Why aren't singularities eliminated for polar plots?
+        elseif (isplot)
+          Z = feval (fun, X);
+          ## Eliminate the singularities
+          Z = __eliminate_sing__ (Z);
+          domain = find_valid_domain (X, [], Z);
+        elseif (iscontour)
+          Z = feval (fun, X, Y);
+          Z = __eliminate_sing__ (Z);
+        else  #  mesh, surf plots
+          Z = feval (fun, X, Y);
+          Z = __eliminate_sing__ (Z);
+          if (circ)
+            ## Use domain calculated at the start.
+            ## The X, Y grids are non-monotonic after conversion from polar
+            ## coordinates and find_valid_domain fails.
+
+          elseif (auto_domain && ! auto_domain_done)
+            valid_domain = find_valid_domain (X, Y, Z);
+            domain_ok = isequal (domain, valid_domain);
+            domain = valid_domain;
+            auto_domain_done = true;  # ensures only 1 round of do loop done
+          else
+            if (! auto_domain_done)
+              domain = find_valid_domain (X, Y, Z);
+            endif
+          end
+        endif
       endif
     endif
-  endif
+  until (domain_ok)
 
+  ## Now, actually call the correct plot function with valid data and domain.
   oldax = gca ();
   unwind_protect
     axes (ax);
     if (iscontour)
-      [clev, h] = feval (pfunc, X, Y, Z);
+      [~, h] = feval (pltfunc, X, Y, Z);
+      ## FIXME: Work around contour setting axis tight.
+      ##        Fix should really be in __countour__.
+      axis (domain);
     elseif (isplot && nargs == 2)
-      h = [];
+      h = zeros (length (XX), 1);
       hold_state = get (ax, "nextplot");
       for i = 1 : length (XX)
-        h = [h; plot(XX{i}, YY{i})];
+        h(i) = plot(XX{i}, YY{i});
         if (i == 1)
           set (ax, "nextplot", "add");
         endif
       endfor
       set (ax, "nextplot", hold_state);
-    elseif (ispolar || isplot)
-      h = feval (pfunc, X, Z);
-      if (isplot && !parametric)
-        axis ([X(1), X(end), yrange]);
+      axis (domain);
+    elseif (isplot || ispolar)
+      h = feval (pltfunc, X, Z);
+      if (isplot && ! parametric)
+        axis (domain);
+      endif
+    elseif (isplot3)
+      if (animate)
+        comet3 (X, Y, Z, .05);  # draw animation, then replace with true plot3
       endif
-    else
-      h = feval (pfunc, X, Y, Z);
+      h = feval (pltfunc, X, Y, Z);
+      set (gca, "box", "off");
+      grid ("on");
+      zlabel ("z");
+    else  # mesh and surf plots
+      h = feval (pltfunc, X, Y, Z);
+      ## FIXME: surf, mesh should really do a better job of setting zlim
+      if (! parametric)
+        axis (domain);
+      endif
     endif
     xlabel (xarg);
     ylabel (yarg);
@@ -439,7 +475,88 @@
 
 endfunction
 
+## Eliminate bad data (complex values, infinities, singularities)
 function x = __eliminate_sing__ (x)
-  x (isinf (x)) = NaN;
-  x (abs (del2 (x)) > 0.2 * (max (x(:)) - min (x(:)))) = NaN;
+  if (iscomplex (x))
+    x(imag (x) != 0) = NaN;
+  endif
+  x(isinf (x)) = NaN;
+  ## High rates of curvature are treated as singularities
+  threshold = 0.2 * (max (x(:)) - min (x(:)));
+  x(abs (del2 (x)) > threshold) = NaN;
 endfunction
+
+## Find: 1) range of function where there are not NaN values,
+##       2) function is changing (not just flat surface)
+function domain = find_valid_domain (X, Y, Z);
+
+  if (isvector (Z))
+    ## 2-D data for isplot
+    domain = [X(1) X(end)];
+
+    ## Guess a range which includes the "mass" of the data by using a 
+    ## median-based approach.  The center 3/4 of the data is used to
+    ## determine the range of the data.
+    ## This seems to be vaguely what Matlab does, but can't be sure.
+    XX = sort (Z(isfinite (Z)));
+    if (length (X) > 4)
+      irlo = XX(fix (1/8 * length (XX)));
+      irhi = XX(fix (7/8 * length (XX)));
+      d = irhi - irlo;
+      domain(3) = max (XX(1) - d/8, irlo - d);
+      domain(4) = min (XX(end) + d/8, irhi + d);
+    else
+      domain(3:4) = [XX(1), XX(end)];
+    endif
+
+    #{
+    ## FIXME: Old algorithm for removing singularities
+    ## Deprecated in 3.8.  Can be removed if no problems appear in ezplot.
+    idx = 2 : length (Z);
+    idx = find (((Z(idx) > yrange(2) / 2) & (Z(idx-1) < yrange(1) / 2)) |
+                ((Z(idx) < yrange(1) / 2) & (Z(idx-1) > yrange(2) / 2)));
+    Z(idx) = NaN;
+    #}
+
+  else
+    ## 3-D data such as mesh, surf
+    Zfinite = ! isnan (Z); 
+    Zrows = any (Zfinite, 2); 
+    rmin = find (Zrows, 1, "first"); 
+    rmax = find (Zrows, 1, "last"); 
+    Zcols = any (Zfinite, 1); 
+    cmin = find (Zcols, 1, "first"); 
+    cmax = find (Zcols, 1, "last"); 
+
+    ## Handle nasty case of all NaNs 
+    if (isempty (rmin))
+      rmin = 1, rmax = rows (Z);
+    endif
+    if (isempty (cmin))
+      cmin = 1, cmax = columns (Z);
+    endif
+
+    if (   ! any (isnan (Z([rmin, rmax],:)(:)))
+        && ! any (isnan (Z(:, [cmin, cmax])(:))))
+      ## Exclude surfaces along borders which are flat (gradient =~ 0).
+      ## Technically, this calculation might be better done with actual
+      ## deltaX, deltaY values.  But, data is usually meshgridded
+      ## (constant spacing) so working with deltaROW#, deltaCOL# is fine.
+      [Zx, Zy] = gradient (Z(rmin:rmax, cmin:cmax));
+      Zgrad = sqrt (Zx.^2 + Zy.^2);
+      slope = ((max (Z(:)) - min (Z(:)))
+                / sqrt ((rmax - rmin)^2 + (cmax - cmin)^2));
+      slope /= 125;  # threshold for discarding points.
+      Zrows = any (Zgrad > slope, 2); 
+      rmin += find (Zrows, 1, "first") - 1; 
+      rmax += find (Zrows, 1, "last") - rows (Zrows); 
+      Zcols = any (Zgrad > slope, 1); 
+      cmin += find (Zcols, 1, "first") - 1; 
+      cmax += find (Zcols, 1, "last") - columns (Zcols); 
+    endif
+
+    domain = [X(1,cmin) X(1,cmax) Y(rmin,1) Y(rmax,1)];
+  endif
+  
+endfunction
+