--- a/scripts/plot/private/__ezplot__.m
+++ b/scripts/plot/private/__ezplot__.m
@@ -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
+