--- a/scripts/plot/fplot.m	Thu Jul 11 13:11:55 2013 +0100
+++ b/scripts/plot/fplot.m	Thu Jul 11 09:25:54 2013 -0700
@@ -20,41 +20,51 @@
 ## @deftypefn  {Function File} {} fplot (@var{fn}, @var{limits})
 ## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{tol})
 ## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{n})
-## @deftypefnx {Function File} {} fplot (@dots{}, @var{fmt})
-## Plot a function @var{fn} within defined limits.
+## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{fmt})
+## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{tol}, @var{n}, @var{fmt})
+## @deftypefnx {Function File} {[@var{x}, @var{y}] =} fplot (@dots{})
+## Plot a function @var{fn} within the range defined by @var{limits}.
+##
 ## @var{fn} is a function handle, inline function, or string
 ## containing the name of the function to evaluate.
-## The limits of the plot are given by @var{limits} of the form
-## @code{[@var{xlo}, @var{xhi}]} or @code{[@var{xlo}, @var{xhi},
-## @var{ylo}, @var{yhi}]}.  @var{tol} is the default tolerance to use for the
-## plot, and if @var{tol} is an integer it is assumed that it defines the
-## number points to use in the plot.  The @var{fmt} argument is passed
-## to the plot command.
+## The limits of the plot are of the form @code{[@var{xlo}, @var{xhi}]} or
+## @code{[@var{xlo}, @var{xhi}, @var{ylo}, @var{yhi}]}.
+## The next three arguments are all optional and any number of them may be
+## given in any order.
+## @var{tol} is the relative tolerance to use for the plot and defaults
+## to 2e-3 (.2%).
+## @var{n} is the minimum number of of points to use.  When @var{n} is
+## specified, the maximum stepsize will be
+## @code{@var{xhi} - @var{xlo} / @var{n}}.  More than @var{n} points may still
+## be used in order to meet the relative tolerance requirement.
+## The @var{fmt} argument specifies the linestyle to be used by the plot
+## command.
+##
+## With no output arguments the results are immediately plotted.  With two
+## output arguments the 2-D plot data is returned.  The data can subsequently
+## be plotted manually with @code{plot (@var{x}, @var{y}).
+##
+## Example:
 ##
 ## @example
 ## @group
-## fplot ("cos", [0, 2*pi])
+## fplot (@cos, [0, 2*pi])
 ## fplot ("[cos(x), sin(x)]", [0, 2*pi])
 ## @end group
 ## @end example
+##
+## Note: @code{fplot} works best with continuous functions.  Functions with
+## discontinuities are unlikely to plot well.  This restriction may be
+## removed in the future.
 ## @seealso{plot}
 ## @end deftypefn
 
 ## Author: Paul Kienzle <pkienzle@users.sf.net>
 
-function fplot (fn, limits, n = 0.002, fmt = "")
-
-  if (nargin < 2 || nargin > 4)
-    print_usage ();
-  endif
+function [X, Y] = fplot (fn, limits, varargin)
 
-  if (iscomplex (limits) || (numel (limits) != 2 && numel (limits) != 4))
-    error ("fplot: LIMITS must be a real vector with 2 or 4 elements");
-  endif
-
-  if (ischar (n))
-    fmt = n;
-    n = 0.002;
+  if (nargin < 2 || nargin > 5)
+    print_usage ();
   endif
 
   if (strcmp (typeinfo (fn), "inline function"))
@@ -71,52 +81,89 @@
     error ("fplot: FN must be a function handle, inline function, or string");
   endif
 
-  if (n != fix (n))
-    tol = n;
-    x0 = linspace (limits(1), limits(2), 5)';
+  if (iscomplex (limits) || (numel (limits) != 2 && numel (limits) != 4))
+    error ("fplot: LIMITS must be a real vector with 2 or 4 elements");
+  endif
+
+  n = 5;
+  tol = 2e-3;
+  fmt = "";
+  for i = 1:numel (varargin)
+    arg = varargin{i};
+    if (ischar (arg))
+      fmt = arg;
+    elseif (isnumeric (arg) && isscalar (arg) && arg > 0)
+      if (arg == fix (arg))
+        n = arg;
+      else
+        tol = arg;
+      endif
+    else
+      error ("fplot: bad input in position %d", i+2);
+    endif
+  endfor
+
+  if (n != 5) 
+    ## n was specified 
+    x0 = linspace (limits(1), limits(2), n/2 + 1)';
     y0 = feval (fn, x0);
-    err0 = Inf;
-    n = 8;
     x = linspace (limits(1), limits(2), n)';
     y = feval (fn, x);
-   
-    if (! size_equal (x0, y0))
-      ## FN is a constant value function
-      y0 = repmat (y0, size (x0));
-      y = repmat (y, size (x));
-    endif
-
-    while (n < 2 .^ 20)
-      y00 = interp1 (x0, y0, x, "linear");
-      err = 0.5 * max (abs ((y00 - y) ./ (y00 + y))(:));
-      if (err == err0 || err < tol)
-        break;
-      endif
-      x0 = x;
-      y0 = y;
-      err0 = err;
-      n = 2 * (n - 1) + 1;
-      x = linspace (limits(1), limits(2), n)';
-      y = feval (fn, x);
-    endwhile
   else
+    x0 = linspace (limits(1), limits(2), 5)';
+    y0 = feval (fn, x0);
+    n = 8;
     x = linspace (limits(1), limits(2), n)';
     y = feval (fn, x);
   endif
 
-  plot (x, y, fmt);
-
-  if (length (limits) > 2)
-    axis (limits);
+  if (rows (x0) != rows (y0))
+    ## FN is a constant value function
+    y0 = repmat (y0, size (x0));
+    y = repmat (y, size (x));
   endif
 
-  if (isvector (y))
-    legend (nam);
+  err0 = Inf;
+
+  ## FIXME: This algorithm should really use adaptive scaling as the 
+  ##        the numerical quadrature algorithms do so that extra points are
+  ##        used where they are needed and not spread evenly over the entire
+  ##        x-range.  Try any function with a discontinuity such as
+  ##        fplot (@tan, [-2, 2]) or fplot ("1./x", [-3, 2]) to see the
+  ##        problems with the current solution.
+
+  while (n < 2^18)    # Something is wrong if we need more than 250K points
+    yi = interp1 (x0, y0, x, "linear");
+    ## relative error calculation using average of [yi,y] as reference
+    ## since neither estimate is known a priori to be better than the other.
+    err = 0.5 * max (abs ((yi - y) ./ (yi + y))(:));
+    if (err < tol || abs (err - err0) < tol/2)
+      ## Either relative tolerance has been met OR
+      ## algorithm has stopped making any reasonable progress per iteration.
+      break;
+    endif
+    x0 = x;
+    y0 = y;
+    err0 = err;
+    n = 2 * (n - 1) + 1;
+    x = linspace (limits(1), limits(2), n)';
+    y = feval (fn, x);
+  endwhile
+
+  if (nargout == 2)
+    X = x;
+    Y = y;
   else
-    for i = 1:columns (y)
-      nams{i} = sprintf ("%s(:,%i)", nam, i);
-    endfor
-    legend (nams{:});
+    plot (x, y, fmt);
+    axis (limits);
+    if (isvector (y))
+      legend (nam);
+    else
+      for i = 1:columns (y)
+        nams{i} = sprintf ("%s(:,%i)", nam, i);
+      endfor
+      legend (nams{:});
+    endif
   endif
 
 endfunction
@@ -124,17 +171,30 @@
 
 %!demo
 %! clf;
-%! fplot ('cos', [0, 2*pi]);
+%! fplot (@cos, [0, 2*pi]);
 
 %!demo
 %! clf;
 %! fplot ('[cos(x), sin(x)]', [0, 2*pi]);
 
+%!demo
+%! clf;
+%! ## sinc function
+%! fh = @(x) sin (x*pi) ./ (x*pi);
+%! fplot (fh, [-5, 5]);
+
+%!test
+%! [x, y] = fplot ("[cos(x), sin(x)]", [0, 2*pi]);
+%! assert (columns (y) == 2);
+%! assert (rows (x) == rows (y));
+%! assert (y, [cos(x), sin(x)], -2e-3);
+
 %% Test input validation
 %!error fplot (1)
-%!error fplot (1,2,3,4,5)
+%!error fplot (1,2,3,4,5,6)
+%!error <FN must be a function handle> fplot (1, [0 1])
 %!error <LIMITS must be a real vector> fplot (@cos, [i, 2*i])
 %!error <LIMITS must be a real vector with 2 or 4> fplot (@cos, [1])
 %!error <LIMITS must be a real vector with 2 or 4> fplot (@cos, [1 2 3])
-%!error <FN must be a function handle> fplot (1, [0 1])
+%!error <bad input in position 3> fplot (@cos,[-1,1], {1})