# HG changeset patch
# User Nicholas R. Jankowski <jankowskin@asme.org>
# Date 1508020251 14400
#      Sat Oct 14 18:30:51 2017 -0400
# Node ID 2c820ab8a0fbc7963e7f158fc2c1c063cf60c9f1
# Parent  d72ba91913527c6ecb6b3a22a4da3326a825cb05
Updated integral2.m and integral3.m and added quad2d (bug #52074)


* /scripts/general/quad2d.m: new function file implementing tiled integration



* /scripts/general/integral2.m: updated to avoid dblquad and make use of quad2d

* /scripts/general/integral3.m: updated to avoid double or triplequad and make use of quad2d
* /scripts/general/module.mk: add quad2d.m to the list

* /scripts/help/__unimplemented__.m: removed quad2d from list and BIST
* quad.txi: updated 2D&3D integration section, added quad2d
* NEWS: announce new function quad2d

diff -r d72ba9191352 -r 2c820ab8a0fb NEWS
--- a/NEWS	Wed Oct 04 14:13:49 2017 -0500
+++ b/NEWS	Sat Oct 14 18:30:51 2017 -0400
@@ -78,6 +78,7 @@
       integral2
       integral3
       openvar
+	  quad2d
       repelem
       rticks
       thetaticks
diff -r d72ba9191352 -r 2c820ab8a0fb doc/interpreter/quad.txi
--- a/doc/interpreter/quad.txi	Wed Oct 04 14:13:49 2017 -0500
+++ b/doc/interpreter/quad.txi	Sat Oct 14 18:30:51 2017 -0400
@@ -274,13 +274,13 @@
 @node Functions of Multiple Variables
 @section Functions of Multiple Variables
 
-Octave does not have built-in functions for computing the integral of
-functions of multiple variables directly.  It is possible, however, to
-compute the integral of a function of multiple variables using the
-existing functions for one-dimensional integrals.
+Octave includes several functions for computing the integral of functions of 
+multiple variables. This procedure can generally be performed by creating a 
+function that integrates @math{f} with repect to @math{x}, and then integrates
+that function with respect to @math{y}. 
+This procedure can be performed manually using the following example which
+integrates the function:
 
-To illustrate how the integration can be performed, we will integrate
-the function
 @tex
 $$
   f(x, y) = \sin(\pi x y)\sqrt{x y}
@@ -295,13 +295,10 @@
 @end ifnottex
 for @math{x} and @math{y} between 0 and 1.
 
-The first approach creates a function that integrates @math{f} with
-respect to @math{x}, and then integrates that function with respect to
-@math{y}.  Because @code{quad} is written in Fortran it cannot be called
-recursively.  This means that @code{quad} cannot integrate a function
-that calls @code{quad}, and hence cannot be used to perform the double
-integration.  Any of the other integrators, however, can be used which is
-what the following code demonstrates.
+Using @code{quadgk) in the example below, a double integration can be 
+performed.  (Note that any of the 1D quadrature functions can be used in this 
+fashion except for @code{quad} since it is written in Fortran and cannot be 
+called recursively.)
 
 @example
 @group
@@ -318,9 +315,9 @@
 @end group
 @end example
 
-The above process can be simplified with the @code{dblquad} and
-@code{triplequad} functions for integrals over two and three variables.  For
-example:
+The above process is implemented in the function @code{dblquad} and the 3D 
+equivalent in @code{triplequad} for integrals over two and three variables.  
+For example:
 
 @example
 @group
@@ -333,19 +330,30 @@
 
 @DOCSTRING(triplequad)
 
-@code{integral2} and @code{integral3} are also provided as partial
-compatibility wrappers to @code{dblquad} and @code{triplequad}:
+A separate 2D integration method is implemented in the function @code{quad2d}. 
+This function performs a 'tiled' integration by subdividing the integration 
+domain into rectangular regions and performing separate integrations over 
+those domains. The domains are further subdivided in areas requiring refinement
+to reach the desired numerical accuracy. For certain functions this method can
+be faster than the 2D iteration used in the other functions above.
+
+@DOCSTRING(quad2d)
+
+Additionally, the functions @code{integral2} and @code{integral3} are  provided
+as general 2D and 3D integration functions. They will autoselect between 
+iterated and tiled integration methods, and unlike @code{dblquad} and 
+@code{triplequad} will work with non-rectangular integration domains. 
 
 @DOCSTRING(integral2)
 
 @DOCSTRING(integral3)
 
-The above mentioned approach works, but is fairly slow, and that problem
-increases exponentially with the dimensionality of the integral.  Another
-possible solution is to use Orthogonal Collocation as described in the
-previous section (@pxref{Orthogonal Collocation}).  The integral of a function
-@math{f(x,y)} for @math{x} and @math{y} between 0 and 1 can be approximated
-using @math{n} points by
+The above integrations can be fairly slow, and that problem increases 
+exponentially with the dimensionality of the integral.  Another possible 
+solution for 2D integration is to use Orthogonal Collocation as described in 
+the previous section (@pxref{Orthogonal Collocation}).  The integral of a 
+function @math{f(x,y)} for @math{x} and @math{y} between 0 and 1 can be 
+approximated using @math{n} points by
 @tex
 $$
  \int_0^1 \int_0^1 f(x,y) d x d y \approx \sum_{i=1}^n \sum_{j=1}^n q_i q_j f(r_i, r_j),
diff -r d72ba9191352 -r 2c820ab8a0fb scripts/general/integral2.m
--- a/scripts/general/integral2.m	Wed Oct 04 14:13:49 2017 -0500
+++ b/scripts/general/integral2.m	Sat Oct 14 18:30:51 2017 -0400
@@ -1,221 +1,302 @@
-## Copyright (C) 2017 Nicholas Jankowski
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 3 of the License, or
-## (at your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-## GNU General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn  {} {@var{q} =} integral2 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb})
-## @deftypefnx {} {@var{q} =} integral2 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{prop}, @var{val}, @dots{})
-##
-## Numerically evaluate the two-dimensional integral of @var{f} using adaptive
-## quadrature over the two-dimensional domain defined by @var{xa}, @var{xb},
-## @var{ya}, @var{yb} (scalars may be finite or infinite).
-##
-## @code{integral2} is a wrapper for @code{dblquad} intended to provide
-## @sc{matlab} compatibility.  More control of the numerical integration may
-## be achievable by calling the various quadrature functions directly.
-##
-## @var{f} is a function handle, inline function, or string containing the name
-## of the function to evaluate.  The function @var{f} must be of the form
-## @math{z = f(x,y)} where @var{x} is a vector and @var{y} is a scalar.  It
-## should return a vector of the same length and orientation as @var{x}.
-##
-## Additional optional parameters can be specified using
-## @qcode{"@var{property}", @var{value}} pairs.  Valid properties are:
-##
-## @table @code
-## @item AbsTol
-## Define the absolute error tolerance for the quadrature.  The default
-## value is 1e-10 (1e-5 for single).
-##
-## @item RelTol
-## Define the relative error tolerance for the quadrature.  The default
-## value is 1e-6 (1e-4 for single).
-## @end table
-##
-## Adaptive quadrature is used to minimize the estimate of error until the
-## following is satisfied:
-## @tex
-## $$error \leq \max \left( AbsTol, RelTol\cdot\vert q\vert \right)$$
-## @end tex
-## @ifnottex
-##
-## @example
-## @group
-##   @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|).
-## @end group
-## @end example
-##
-## @end ifnottex
-##
-## Known @sc{matlab} incompatibilities:
-##
-## @enumerate
-## @item
-## If tolerances are left unspecified, and any integration limits or waypoints
-## are of type @code{single}, then Octave's integral functions automatically
-## reduce the default absolute and relative error tolerances as specified
-## above.  If tighter tolerances are desired they must be specified.
-## @sc{matlab} leaves the tighter tolerances appropriate for @code{double}
-## inputs in place regardless of the class of the integration limits.
-#
-## @item
-## @code{integral2} currently only integrates functions over rectangular
-## domains.  Implementing @var{ya} and @var{yb} as functions of @var{x} is a
-## planned future improvement.
-##
-## @item
-## The @qcode{"Method"} property is not yet implemented in Octave due to the
-## lack of a @qcode{"tiled"} integrator implementation.  All integrals are
-## evaluated using an equivalent of the @qcode{"iterated"} method.
-## @end enumerate
-##
-## @seealso{integral, integral3, quad, quadgk, quadv, quadl, quadcc, trapz,
-##          dblquad, triplequad}
-## @end deftypefn
-
-function q = integral2 (f, xa, xb, ya, yb, varargin)
-  ## FIXME: It is possible that a non-rectangular domain could be handled by
-  ##        overlaying the integrand with a boolean mask function such that
-  ##        the integration occurs over a rectangle, but regions outside the
-  ##        desired domain contribute zero to the integral.  This may be an
-  ##        inefficient but acceptable hack to get around the rectangular
-  ##        domain limit without having to rewrite the integrating function.
-
-  ## FIXME: Implement "Method" property to let the user select between iterated
-  ##        and tiled integration.  Tiled integration follows the method of
-  ##        Matlab's quad2d function, currently unimplemented in Octave.
-  ##        Should probably just wait for a quad2d implementation to point the
-  ##        integral2 wrapper to, instead of trying to re-create it here.  The
-  ##        following can be added to the help docstring once it is functional:
-  ## @item Method
-  ## Specifies the two dimensional integration method to be used, with valid
-  ## options being @var{"auto"}, @var{"tiled"}, or @var{"iterated"}.
-  ## @code{integral} will use @var{"auto"} by default, where it will usually
-  ## choose @var{"tiled"} unless any of the integration limits are infinite.
-
-  if (nargin < 5 || (mod (nargin, 2) == 0))
-    print_usage ();
-  endif
-
-  if (! is_function_handle (f))
-    print_usage ();
-  endif
-
-  if (! (isscalar (xa) && isscalar (xb) && isscalar (ya)) && isscalar (yb))
-    print_usage ();
-  endif
-
-  ## Check for single or double limits to set appropriate default tolerance.
-  issingle = isa ([xa, xb, ya, yb], "single");
-
-  ## Set defaults, update with any specified parameters.
-  if (issingle)
-    abstol = 1e-5;
-    reltol = 1e-4;
-  else
-    abstol = 1e-10;
-    reltol = 1e-6;
-  endif
-
-  if (nargin == 5)
-    ## Pass the simplest case directly to integrator.
-    ## Let quadcc function handle input checks.
-    q = dblquad (f, xa, xb, ya, yb, [abstol, reltol], @quadcc);
-
-  else
-    ## Parse options to determine how to call integrator.
-    intmethod = [];
-
-    idx = 1;
-    while (idx < nargin - 5)
-      prop = varargin{idx++};
-      if (! ischar (prop))
-        error ("integral2: property PROP must be a string");
-      endif
-
-      switch (tolower (prop))
-        case "abstol"
-          abstol = varargin{idx++};
-          if (! (isnumeric (abstol) && isscalar (abstol) && abstol >= 0))
-            error ("integral2: AbsTol value must be a numeric scalar >= 0");
-          endif
-
-        case "reltol"
-          reltol = varargin{idx++};
-          if (! (isnumeric (reltol) && isscalar (reltol) && reltol >= 0))
-            error ("integral2: RelTol value must be a numeric scalar >= 0");
-          endif
-
-        case "method"
-          intmethod = varargin{idx++};
-          warning (["integral2: Only 'iterated' method implemented.  ", ...
-                    "Method property ignored."]);
-        otherwise
-          error ("integral2: unknown property '%s'", prop);
-      endswitch
-    endwhile
-
-    q = dblquad (f, xa, xb, ya, yb, [abstol, reltol], @quadcc);
-
-  endif
-
-endfunction
-
-
-%!test
-%! f = @(x, y) x .* y;
-%! assert (integral2 (f, 0, 1, 0, 1), 0.25, 1e-10);
-
-%!test
-%! f = @(x, y) 9 * x.^2 + 15 * y.^2;
-%! assert (integral2 (f, 0, 5, -5, 0, "AbsTol", 1e-9), 5000, 1e-9);
-%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6), 5000, -1e-6);
-%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6, "AbsTol", 1e-9),
-%!         5000, 1e-9);
-
-## tests from dblquad
-%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-7),
-%!        2*log (2), 1e-7)
-%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "RelTol", 1e-6),
-%!        2*log (2), -1e-6)
-%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-8,
-%!        "RelTol", 1e-6), 2*log (2), -1e-6)
-%!assert (integral2 (@(x, y) exp (-x.^2 - y.^2) , -1, 1, -1, 1),
-%!        pi * erf (1).^2, 1e-10)
-
-%!assert (integral2 (@plus, 1, 2, 3, 4), 5, 1e-10)
-
-## Test input validation
-%!error integral2
-%!error integral2 (0, 1 ,2 ,3 ,4)
-%!error integral2 (@plus)
-%!error integral2 (@plus, 1)
-%!error integral2 (@plus, 1, 2)
-%!error integral2 (@plus, 1, 2, 3)
-%!error integral2 (@plus, 1, 2, 3, [4 5])
-%!error integral2 (@plus, 1, 2, 3, "test")
-%!error integral2 (@plus, 1, 2, 3, 4, "foo")
-%!error <unknown property 'foo'>  integral2 (@plus, 1, 2, 3, 4, "foo", "bar")
-%!error <property PROP must be a string> integral2 (@plus, 1, 2, 3, 4, 99, "bar")
-%!error <AbsTol value must be a numeric> integral2 (@plus, 1, 2, 3, 4, "AbsTol", "foo")
-%!error <AbsTol value must be a .* scalar> integral2 (@plus, 1, 2, 3, 4, "AbsTol", [1, 2])
-%!error <AbsTol value must be.* .= 0> integral2 (@plus, 1, 2, 3, 4, "AbsTol", -1)
-%!error <RelTol value must be a numeric> integral2 (@plus, 1, 2, 3, 4, "RelTol", "foo")
-%!error <RelTol value must be a .* scalar> integral2 (@plus, 1, 2, 3, 4, "RelTol", [1, 2])
-%!error <RelTol value must be.* .= 0> integral2 (@plus, 1, 2, 3, 4, "RelTol", -1)
-%!warning <Only 'iterated' method implemented>
-%! q = integral2 (@plus, 0, 1, 0, 1, "Method", "tiled");
+## Copyright (C) 2017 Nicholas Jankowski, David Bateman
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {} {@var{q} =} integral2 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb})
+## @deftypefnx {} {@var{q} =} integral2 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{prop}, @var{val}, @dots{})
+## @deftypefnx {} {[@var{q}, @var{err}] =} integral2 (@dots{})
+##
+## Numerically evaluate the two-dimensional integral of @var{f} using adaptive
+## quadrature over the two-dimensional domain defined by @var{xa}, @var{xb},
+## @var{ya}, @var{yb} (scalars may be finite or infinite). Additionally, 
+## @var{ya} and @var{yb} may be scalar functions of @var{x}, allowing for the
+## integration over non-rectangular domains.
+##
+## @var{f} is a function handle, inline function, or string containing the name
+## of the function to evaluate.  The function @var{f} must be of the form
+## @math{z = f(x,y)} where @var{x} is a vector and @var{y} is a scalar.  It
+## should return a vector of the same length and orientation as @var{x}.
+##
+## Additional optional parameters can be specified using
+## @qcode{"@var{property}", @var{value}} pairs.  Valid properties are:
+##
+## @table @code
+## @item Method
+## Specifies the two dimensional integration method to be used, with valid
+## options being @var{"auto"}, @var{"tiled"}, or @var{"iterated"}.
+## @code{integral} will use @var{"auto"} by default, where it will usually
+## choose @var{"tiled"} unless any of the integration limits are infinite.
+##
+## @item Vectorized
+## Option to enable/disable vectorized integration. False forces octave to use 
+## only scalar inputs when calling the integrand, which enables integrands 
+## @math{f(x,y)} that have not been vectorized and only accept @var{x} and 
+## @var{y} as scalars. Default value is @code{true}.
+##
+## @item AbsTol
+## Define the absolute error tolerance for the quadrature.  The default
+## value is 1e-10 (1e-5 for single).
+##
+## @item RelTol
+## Define the relative error tolerance for the quadrature.  The default
+## value is 1e-6 (1e-4 for single).
+## @end table
+##
+## Adaptive quadrature is used to minimize the estimate of error until the
+## following is satisfied:
+## @tex
+## $$error \leq \max \left( AbsTol, RelTol\cdot\vert q\vert \right)$$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+##   @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|).
+## @end group
+## @end example
+##
+## @end ifnottex
+##
+## @var{err} is an approximate bound on the error in the integral
+## @code{abs (@var{q} - @var{I})}, where @var{I} is the exact value of the
+## integral.
+##
+## Known @sc{matlab} incompatibilities:
+##
+## @enumerate
+## @item
+## If tolerances are left unspecified, and any integration limits
+## are of type @code{single}, then Octave's integral functions automatically
+## reduce the default absolute and relative error tolerances as specified
+## above.  If tighter tolerances are desired they must be specified.
+## @sc{matlab} leaves the tighter tolerances appropriate for @code{double}
+## inputs in place regardless of the class of the integration limits.
+## @end enumerate
+##
+##
+## Reference: @nospell{L.F. Shampine},
+## @cite{"@sc{matlab} program for quadrature in 2D"}, Applied Mathematics and 
+## Computation, pp. 266--274, Vol 1, 2008.
+##
+## @seealso{integral, integral3, quad, quadgk, quadv, quadl, quadcc, trapz,
+##          quad2d, dblquad, triplequad}
+## @end deftypefn
+
+function [q, err] = integral2 (f, xa, xb, ya, yb, varargin)
+  
+  if ((nargin < 5) || (mod (nargin, 2) == 0))
+    print_usage ();
+  endif
+
+  if (! is_function_handle (f))
+    print_usage ();
+  endif
+
+  if (! (isscalar (xa) && isscalar (xb)))
+    print_usage ();
+  endif
+
+  ## Check for single or double limits to set appropriate default tolerance.
+  issingle = isa ([xa, xb], "single") || ...
+                   ((! is_function_handle (ya)) && isa (ya, "single")) || ...
+                   ((! is_function_handle (yb)) && isa (yb, "single"));
+  
+  
+  ## Set defaults, update with any specified parameters.
+  if issingle
+    abstol = 1e-5;
+    reltol = 1e-4;
+  else
+    abstol = 1e-10;
+    reltol = 1e-6;
+  endif
+  
+  method = "auto";
+  idx = 1;
+  while (idx < nargin - 5)
+    prop = varargin{idx++};
+    if (! ischar (prop))
+      error ("integral2: property PROP must be a string");
+    endif
+
+    switch (tolower (prop))
+      case "abstol"
+        abstol = varargin{idx++};
+        if (! ((isnumeric (abstol)) && (isscalar (abstol)) && (abstol >= 0)))
+          error ("integral2: AbsTol value must be a numeric scalar >= 0");
+        endif
+
+      case "reltol"
+        reltol = varargin{idx++};
+        if (! ((isnumeric (reltol)) && (isscalar (reltol)) && (reltol >= 0)))
+          error ("integral2: RelTol value must be a numeric scalar >= 0");
+        endif
+
+      case "method"
+        method = tolower (varargin{idx++});
+        if (! any (strcmp (method, {"auto", "iterated", "tiled"})))
+          error ("integral2 : method '%s' unrecognized", method);
+        endif
+      case "vectorized"
+        # option to allow unvectorized functions to be used
+        vectorized = varargin{idx++};
+        if (! islogical (vectorized))
+          error ("integral2: 'vectorized' must be a logical value");
+        elseif (! vectorized)
+		  f = @(x, y) arrayfun (f, x, y);
+		endif      
+      
+      otherwise
+        error ("integral2: unknown property '%s'", prop);
+    endswitch
+  endwhile
+    
+  if strcmp (method, "auto")
+    if ((isinf (xa)) || (isinf (xb)) || ...
+         ((! is_function_handle (ya)) && (isinf (ya))) || ...
+         ((! is_function_handle (yb)) && (isinf (yb))))
+      method = "iterated";
+    else
+      method = "tiled";
+    endif
+  endif
+
+  # check upper and lower bounds of y
+  if (! is_function_handle (ya))
+    if (isscalar (ya))
+      ya = @(x) ya * ones (rows (x), columns (x));
+    else
+      error ("integral2: 'ya' must be a constant or a (vectorized) function.");
+    endif
+  endif
+  if (! is_function_handle (yb))
+    if (isscalar (yb))
+      yb = @(x) yb * ones (rows (x), columns (x));
+    else
+      error ("integral2: 'ya' must be a constant or a (vectorized) function.");
+    endif
+  endif
+  
+  if strcmp (method, "iterated") 
+    q = outer_iterated (f, xa, xb, ya, yb, abstol, reltol);
+    
+    if nargout == 2
+      warning ("integral2: 'iterated' method can not return estimated error");
+      err = 0;
+    endif
+    
+  else
+    [q, err] = quad2d (f, xa, xb, ya, yb, 'AbsTol', abstol, 'RelTol', reltol);
+  endif  
+endfunction
+
+function q = outer_iterated (f, xa, xb, ya, yb, abstol, reltol)   
+    inner = @inner_iterated;
+    q = feval (@quadcc, @(x) inner (x, f, ya, yb, abstol, reltol), xa, xb, ...
+                [abstol, reltol]);
+endfunction
+
+function q = inner_iterated (x, f, ya, yb, abstol, reltol)
+    q = zeros (size (x));
+    for i = 1 : length (x)
+      q(i) = feval (@quadcc, @(y) f(x(i), y), ya(x(i)), yb(x(i)), ...
+                     [abstol, reltol]);
+    endfor
+endfunction
+
+## method tests
+%!test
+%! f = @(x, y) x .* y;
+%! assert (integral2 (f, 0, 1, 0, 1), 0.25, 1e-10);
+%! assert (integral2 (f, 0, 1, 0, 1, "method", "tiled"), 0.25, 1e-10);
+%! assert (integral2 (f, 0, 1, 0, 1, "method", "iterated"), 0.25, 1e-10);
+%! assert (integral2 (f, 0, 1, 0, 1, "method", "auto"), 0.25, 1e-10);
+
+## vectorized = false test
+%!test
+%! f = @(x, y) x * y;
+%! assert (integral2 (f, 0, 1, 0, 1, "vectorized", false), 0.25, 1e-10);
+
+## tolerance tests
+%!test
+%! f = @(x, y) 9 * x.^2 + 15 * y.^2;
+%! assert (integral2 (f, 0, 5, -5, 0, "AbsTol", 1e-9), 5000, 1e-9);
+%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6), 5000, -1e-6);
+%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6, "AbsTol", 1e-9),
+%!         5000, 1e-9);
+
+## tests from dblquad
+%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-7),
+%!        2*log (2), 1e-6);
+%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "RelTol", 1e-6),
+%!        2*log (2), -1e-6);
+%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-8,
+%!        "RelTol", 1e-6), 2*log (2), -1e-6);
+%!assert (integral2 (@(x, y) exp (-x.^2 - y.^2) , -1, 1, -1, 1),
+%!        pi * erf (1).^2, 1e-10);
+
+%!assert (integral2 (@plus, 1, 2, 3, 4), 5, 1e-10);
+
+%!assert (integral2 (@(x,y) 1 ./ (x + y), 0, 1, 0, @(x) 1 - x), 1, -1e-6);
+
+% method="iterated" 
+%!test
+%! f = @(x, y) x .* y;
+%! assert (integral2 (f, 0, 1, 0, 1, "method", "iterated"), 0.25, 1e-10);
+
+%!test
+%! f = @(x, y) 9 * x.^2 + 15 * y.^2;
+%! assert (integral2 (f, 0, 5, -5, 0, "AbsTol", 1e-9, "method", "iterated"), 
+%!         5000, 1e-9);
+%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6, "method", "iterated"), 
+%!         5000, -1e-6);
+%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6, "method", "iterated",
+%!         "AbsTol", 1e-9),5000, 1e-9);
+
+## tests from dblquad
+%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-7, "method",
+%!         "iterated"), 2*log (2), 1e-7);
+%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "RelTol", 1e-6, "method",
+%!         "iterated"), 2*log (2), -1e-6);
+%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-8,
+%!        "RelTol", 1e-6, "method", "iterated"), 2*log (2), -1e-6);
+%!assert (integral2 (@(x, y) exp (-x.^2 - y.^2) , -1, 1, -1, 1, "method",
+%!          "iterated"), pi * erf (1).^2, 1e-10);
+%!assert (integral2 (@plus, 1, 2, 3, 4, "method", "iterated"), 5, 1e-10);
+%!assert (integral2 (@(x,y) 1 ./ (x + y), 0, 1, 0, @(x) 1 - x, "method", 
+%!         "iterated"), 1, -1e-6);
+
+## Test input validation
+%!error integral2
+%!error integral2 (0, 1 ,2 ,3 ,4)
+%!error integral2 (@plus)
+%!error integral2 (@plus, 1)
+%!error integral2 (@plus, 1, 2)
+%!error integral2 (@plus, 1, 2, 3)
+%!error integral2 (@plus, 1, 2, 3, [4 5])
+%!error integral2 (@plus, 1, 2, 3, "test")
+%!error integral2 (@plus, 1, 2, 3, 4, "foo")
+%!error <unknown property 'foo'>  integral2 (@plus, 1, 2, 3, 4, "foo", "bar")
+%!error <property PROP must be a string> integral2 (@plus, 1, 2, 3, 4, 99, "bar")
+%!error <AbsTol value must be a numeric> integral2 (@plus, 1, 2, 3, 4, "AbsTol", "foo")
+%!error <AbsTol value must be a .* scalar> integral2 (@plus, 1, 2, 3, 4, "AbsTol", [1, 2])
+%!error <AbsTol value must be.* .= 0> integral2 (@plus, 1, 2, 3, 4, "AbsTol", -1)
+%!error <RelTol value must be a numeric> integral2 (@plus, 1, 2, 3, 4, "RelTol", "foo")
+%!error <RelTol value must be a .* scalar> integral2 (@plus, 1, 2, 3, 4, "RelTol", [1, 2])
+%!error <RelTol value must be.* .= 0> integral2 (@plus, 1, 2, 3, 4, "RelTol", -1)
+%!error <method 'bad' unrecognized> integral2 (@plus, 1, 2, 3, 4, "method", "bad")
diff -r d72ba9191352 -r 2c820ab8a0fb scripts/general/integral3.m
--- a/scripts/general/integral3.m	Wed Oct 04 14:13:49 2017 -0500
+++ b/scripts/general/integral3.m	Sat Oct 14 18:30:51 2017 -0400
@@ -1,217 +1,290 @@
-## Copyright (C) 2017 Nicholas Jankowski
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 3 of the License, or
-## (at your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-## GNU General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn  {} {@var{q} =} integral3 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb})
-## @deftypefnx {} {@var{q} =} integral3 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}, @var{prop}, @var{val}, @dots{})
-##
-## Numerically evaluate the three-dimensional integral of @var{f} using
-## adaptive quadrature over the three-dimensional domain defined by
-## @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb} (scalars may
-## be finite or infinite).
-##
-## @code{integral3} is a wrapper for @code{triplequad} intended to provide
-## @sc{matlab} compatibility.  More control of the numerical integration may be
-## achievable by calling the various quadrature functions directly.
-##
-## @var{f} is a function handle, inline function, or string containing the name
-## of the function to evaluate.  The function @var{f} must be of the form
-## @math{w = f(x,y,z)} where either @var{x} or @var{y} is a vector and the
-## remaining inputs are scalars.  @var{f} should return a vector of the same
-## length and orientation as the vector input @var{x} or @var{y}.
-##
-## Additional optional parameters can be specified using
-## @qcode{"@var{property}", @var{value}} pairs.  Valid properties are:
-##
-## @table @code
-## @item AbsTol
-## Define the absolute error tolerance for the quadrature.  The default
-## value is 1e-10 (1e-5 for single).
-##
-## @item RelTol
-## Define the relative error tolerance for the quadrature.  The default
-## value is 1e-6 (1e-4 for single).
-## @end table
-##
-## Adaptive quadrature is used to minimize the estimate of error until the
-## following is satisfied:
-## @tex
-## $$error \leq \max \left( AbsTol, RelTol\cdot\vert q\vert \right)$$
-## @end tex
-## @ifnottex
-##
-## @example
-## @group
-##   @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|).
-## @end group
-## @end example
-##
-## @end ifnottex
-##
-## Known @sc{matlab} incompatibilities:
-##
-## @enumerate
-## @item
-## If tolerances are left unspecified, and any integration limits or waypoints
-## are of type @code{single}, then Octave's integral functions automatically
-## reduce the default absolute and relative error tolerances as specified
-## above.  If tighter tolerances are desired they must be specified.
-## @sc{matlab} leaves the tighter tolerances appropriate for @code{double}
-## inputs in place regardless of the class of the integration limits.
-##
-## @item
-## @code{integral3} currently only integrates functions over rectangular
-## volumes.  Implementing @var{ya} and @var{yb} as functions of @var{x}, and
-## @var{za} and @var{zb} as functions of (@var{x,y}) is a planned future
-## improvement.
-##
-## @item
-## The @qcode{"Method"} property is not yet implemented in Octave due to the
-## lack of a @qcode{"tiled"} integrator implementation.  All integrals are
-## evaluated using an equivalent of the @qcode{"iterated"} method.
-## @end enumerate
-##
-## @seealso{integral, integral2, quad, quadgk, quadv, quadl, quadcc, trapz,
-##          dblquad, triplequad}
-## @end deftypefn
-
-function q = integral3 (f, xa, xb, ya, yb, za, zb, varargin)
-  ## FIXME: it is possible that a non-rectangular domain could be handled by
-  ##        overlaying the integrand with a boolean mask function such that
-  ##        the integration occurs over a rectangle but regions outside the
-  ##        desired domain contribute zero to the integral. This may be an
-  ##        inefficient but acceptable hack to get around the rectangular domain
-  ##        limit without having to rewrite the integrating function.
-
-  ## FIXME: implement 'method' property to let the user select between iterated
-  ##        and tiled integration. Tiled integration follows the method of
-  ##        matlab's quad2d function, currently unimplemented in Octave. Should
-  ##        probably just wait for a quad2d implementation to point the
-  ##        integral3 wrapper to instead of trying to recreate it here. The
-  ##        following can be added to the help docstring once it is functional:
-  ## @item Method
-  ## Specifies the underlying 2D integration method to be used on the y and z
-  ## dimensions, with valid options being @var{"auto"}, @var{"tiled"}, or
-  ## @var{"iterated"}. @code{integral3} will use @var{"auto"} by default, where
-  ## it will usually choose @var{"tiled"} unless any of the integration limits
-  ## are infinite.
-
-  if (nargin < 7 || (mod (nargin, 2) == 0))
-    print_usage ();
-  endif
-
-  if (! is_function_handle (f))
-    print_usage ();
-  endif
-
-  if (! (isscalar (xa) && isscalar (xb)
-         && isscalar (ya) && isscalar (yb)
-         && isscalar (za) && isscalar (zb)))
-    print_usage ();
-  endif
-
-  ## Check for single or double limits to set appropriate default tolerance.
-  issingle = isa ([xa, xb, ya, yb, za, zb], "single");
-
-  ## Set defaults, update with any specified parameters.
-  if (issingle)
-    abstol = 1e-5;
-    reltol = 1e-4;
-  else
-    abstol = 1e-10;
-    reltol = 1e-6;
-  endif
-
-  if (nargin == 7)
-    ## Pass the simplest case directly to integrator.
-    ## Let quadcc function handle input checks.
-    q = triplequad (f, xa, xb, ya, yb, za, zb, [abstol, reltol], @quadcc);
-
-  else
-    ## Parse options to determine how to call integrator.
-    intmethod = [];
-
-    idx = 1;
-    while (idx < nargin - 7)
-      prop = varargin{idx++};
-      if (! ischar (prop))
-        error ("integral3: property PROP must be a string");
-      endif
-
-      switch (tolower (prop))
-        case "abstol"
-          abstol = varargin{idx++};
-          if (! (isnumeric (abstol) && isscalar (abstol) && abstol >= 0))
-            error ("integral3: AbsTol value must be a numeric scalar >= 0");
-          endif
-
-        case "reltol"
-          reltol = varargin{idx++};
-          if (! (isnumeric (reltol) && isscalar (reltol) && reltol >= 0))
-            error ("integral2: RelTol value must be a numeric scalar >= 0");
-          endif
-
-        case "method"
-          intmethod = varargin{idx++};
-          warning (["integral3: Only 'iterated' method implemented.  ", ...
-                    "Method property ignored."]);
-        otherwise
-          error ("integral3: unknown property '%s'", prop);
-      endswitch
-    endwhile
-
-    q = triplequad (f, xa, xb, ya, yb, za, zb, [abstol, reltol], @quadcc);
-
-  endif
-
-endfunction
-
-
-%!test
-%! f = @(x, y, z) x.*y.*z;
-%! assert (integral3 (f, 0, 1, 0, 1, 0, 1), 0.125, 1e-10);
-
-%!test
-%! f = @(x,y,z) y.*sin(x) + z.*cos(x);
-%! assert (integral3 (f, 0, pi, 0, 1, -1, 1), 2, 1e-10);
-
-## tests from triplequad
-%! assert (integral3 (@(x,y,z) exp (-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1,
-%!         1), pi^(3/2) * erf (1).^3, 1e-10);
-
-## Test input validation
-%!error integral3 (0, 1 ,2 ,3 ,4, 5, 6)
-%!error integral3 (@plus)
-%!error integral3 (@plus, 1)
-%!error integral3 (@plus, 1, 2)
-%!error integral3 (@plus, 1, 2, 3)
-%!error integral3 (@plus, 1, 2, 3, 4)
-%!error integral3 (@plus, 1, 2, 3, 4, 5)
-%!error integral3 (@plus, 1, 2, 3, 4, 5, [6 7])
-%!error integral3 (@plus, 1, 2, 3, 4, 5, "test")
-%!error integral3 (@plus, 1, 2, 3, 4, 5, 6, "foo")
-%!error <unknown property 'foo'> integral3 (@plus, 1, 2, 3, 4, 5, 6, "foo", "bar")
-%!error <property PROP must be a string> integral3 (@plus, 1, 2, 3, 4, 5, 6, NA, "bar")
-%!error <AbsTol value must be a numeric> integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", "foo")
-%!error <AbsTol value must be a .* scalar> integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", [1, 2])
-%!error <AbsTol value must be.* .= 0> integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", -1)
-%!error <RelTol value must be a numeric> integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", "foo")
-%!error <RelTol value must be a .* scalar> integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", [1, 2])
-%!error <RelTol value must be.* .= 0> integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", -1)
-%!warning <Only 'iterated' method implemented>
-%! q = integral3 (@plus, 0, 1, 0, 1, 0, 1, "Method", "tiled");
+## Copyright (C) 2017 David Bateman
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {} {@var{q} =} integral3 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb})
+## @deftypefnx {} {@var{q} =} integral3 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}, @var{prop}, @var{val}, @dots{})
+##
+## Numerically evaluate the three-dimensional integral of @var{f} using adaptive
+## quadrature over the three-dimensional domain defined by @var{xa}, @var{xb},
+## @var{ya}, @var{yb}, @var{za}, @var{zb} (scalars may be finite or infinite). Additionally, 
+## @var{ya} and @var{yb} may be scalar functions of @var{x} and @var{za} and @var{zb}
+## maybe be scalar functions of @var{x} and @var{y}, allowing for the
+## integration over non-rectangular domains.
+##
+## @var{f} is a function handle, inline function, or string containing the name
+## of the function to evaluate.  The function @var{f} must be of the form
+## @math{z = f(x,y)} where @var{x} is a vector and @var{y} is a scalar.  It
+## should return a vector of the same length and orientation as @var{x}.
+##
+## Additional optional parameters can be specified using
+## @qcode{"@var{property}", @var{value}} pairs.  Valid properties are:
+##
+## @table @code
+## @item Method
+## Specifies the two dimensional integration method to be used, with valid
+## options being @var{"auto"}, @var{"tiled"}, or @var{"iterated"}.
+## @code{integral} will use @var{"auto"} by default, where it will usually
+## choose @var{"tiled"} unless any of the integration limits are infinite.
+##
+## @item Vectorized
+## Option to enable/disable vectorized integration. False forces octave to use 
+## only scalar inputs when calling the integrand, which enables integrands 
+## @math{f(x,y)} that have not been vectorized and only accept @var{x} and 
+## @var{y} as scalars. Default value is @code{true}.
+##
+## @item AbsTol
+## Define the absolute error tolerance for the quadrature.  The default
+## value is 1e-10 (1e-5 for single).
+##
+## @item RelTol
+## Define the relative error tolerance for the quadrature.  The default
+## value is 1e-6 (1e-4 for single).
+## @end table
+##
+## Adaptive quadrature is used to minimize the estimate of error until the
+## following is satisfied:
+## @tex
+## $$error \leq \max \left( AbsTol, RelTol\cdot\vert q\vert \right)$$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+##   @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|).
+## @end group
+## @end example
+##
+## @end ifnottex
+##
+## @var{err} is an approximate bound on the error in the integral
+## @code{abs (@var{q} - @var{I})}, where @var{I} is the exact value of the
+## integral.
+##
+## Known @sc{matlab} incompatibilities:
+##
+## @enumerate
+## @item
+## If tolerances are left unspecified, and any integration limits
+## are of type @code{single}, then Octave's integral functions automatically
+## reduce the default absolute and relative error tolerances as specified
+## above.  If tighter tolerances are desired they must be specified.
+## @sc{matlab} leaves the tighter tolerances appropriate for @code{double}
+## inputs in place regardless of the class of the integration limits.
+## @end enumerate
+##
+##
+## Reference: @nospell{L.F. Shampine},
+## @cite{"@sc{matlab} program for quadrature in 2D"}, Applied Mathematics and 
+## Computation, pp. 266--274, Vol 1, 2008.
+##
+## @seealso{integral, integral2, quad, quadgk, quadv, quadl, quadcc, trapz,
+##          quad2d, dblquad, triplequad}
+## @end deftypefn
+
+function q = integral3 (f, xa, xb, ya, yb, za, zb, varargin)
+  
+  if (nargin < 7 || (mod (nargin, 2) == 0))
+    print_usage ();
+  endif
+
+  if (! is_function_handle (f))
+    print_usage ();
+  endif
+
+  if (! (isscalar (xa) && isscalar (xb)))
+    print_usage ();
+  endif
+
+  ## Check for single or double limits to set appropriate default tolerance.
+  issingle = isa ([xa, xb], "single") || ...
+             ((! is_function_handle (ya)) && isa (ya, "single")) || ...
+             ((! is_function_handle (yb)) && isa (yb, "single")) || ...
+             ((! is_function_handle (za)) && isa (za, "single")) || ...
+             ((! is_function_handle (zb)) && isa (zb, "single"));
+  
+  ## Set defaults, update with any specified parameters.
+  if issingle
+    abstol = 1e-5;
+    reltol = 1e-4;
+  else
+    abstol = 1e-10;
+    reltol = 1e-6;
+  endif
+  
+  method = "auto";
+  vectorized = true;
+  idx = 1;
+  while (idx < nargin - 7)
+    prop = varargin{idx++};
+    if (! ischar (prop))
+      error ("integral3: property PROP must be a string");
+    endif
+
+    switch (tolower (prop))
+      case "abstol"
+        abstol = varargin{idx++};
+        if (! ((isnumeric (abstol)) && (isscalar (abstol)) && (abstol >= 0)))
+          error ("integral3: AbsTol value must be a numeric scalar >= 0");
+        endif
+
+      case "reltol"
+        reltol = varargin{idx++};
+        if (! ((isnumeric (reltol)) && (isscalar (reltol)) && (reltol >= 0)))
+          error ("integral3: RelTol value must be a numeric scalar >= 0");
+        endif
+
+      case "method"
+        method = tolower (varargin{idx++});
+        if (! any (strcmp (method, {"auto", "iterated", "tiled"})))
+          error ("integral3 : method '%s' unrecognized", method);
+        endif
+      case "vectorized"
+        # option to allow unvectorized functions to be used
+        vectorized = varargin{idx++};
+        if (! islogical (vectorized))
+          error ("integral3: 'vectorized' must be a logical value");
+        endif
+      otherwise
+        error ("integral3: unknown property '%s'", prop);
+    endswitch
+  endwhile
+    
+  if strcmp (method, "auto")
+    if ((isinf (xa)) || (isinf (xb)) || ...
+         ((! is_function_handle(ya)) && (isinf (ya))) || ...
+         ((! is_function_handle(yb)) && (isinf (yb))) || ...
+         ((! is_function_handle(za)) && (isinf (za))) || ...
+         ((! is_function_handle(zb)) && (isinf (zb))))         
+      method = "iterated";
+    else
+      method = "tiled";
+    endif
+  endif
+
+  # check upper and lower bounds of y
+  if (! is_function_handle (ya))
+    if isscalar (ya)
+      ya = @(x) ya * ones (size (x));
+    else
+      error ("integral3: 'ya' must be a constant or a (vectorized) function.");
+    endif
+  endif
+  if (! is_function_handle (yb))
+    if isscalar (yb)
+      yb = @(x) yb * ones (size (x));
+    else
+      error ("integral3: 'ya' must be a constant or a (vectorized) function.");
+    endif
+  endif
+  
+  # check upper and lower bounds of z
+  if (! is_function_handle (za))
+    if isscalar (za)
+      za = @(x, y) za * ones (size(y));
+    else
+      error ("integral3: 'za' must be a constant or a (vectorized) function.");
+    endif
+  endif
+  if (! is_function_handle (zb))
+    if isscalar (zb)
+      zb = @(x, y) zb * ones (size (y));
+    else
+      error ("integral3: 'za' must be a constant or a (vectorized) function.")
+    endif
+  endif
+  
+  inner = @inner;
+  q = feval (@quadcc, @(x) inner (x, f, ya, yb, za, zb, vectorized, method, ...
+          abstol, reltol), xa, xb, [abstol, reltol]);
+
+endfunction
+
+function q = inner (x, f, ya, yb, za, zb, vectorized, method, abstol, reltol)
+    q = zeros (size (x));
+    for i = 1 : length (x)
+      za2 = @(y) za(x(i), y);      
+      zb2 = @(y) zb(x(i), y);
+      f2 = @(y, z) f(x(i), y, z);
+      if (! vectorized)
+        f2 = @(y, z) arrayfun (f2, y, z);
+      endif
+      if strcmp (method, "iterated")
+        inner = @inner_iterated; 
+        q(i) = feval (@quadcc, @(y) inner (y, f2, za2, zb2, abstol, reltol), ...
+              ya(x(i)), yb(x(i)), [abstol, reltol]);      
+      else
+        q(i) = quad2d (f2, ya(x(i)), yb(x(i)), za2, zb2, 'AbsTol', abstol, ...
+                       'RelTol', reltol);
+      endif
+    endfor  
+endfunction
+
+function q = inner_iterated (y, f2, za2, zb2, abstol, reltol)
+    q = zeros (size (y));
+    for i = 1 : length (y)
+      q(i) = feval (@quadcc, @(z) f2(y(i), z), za2(y(i)), zb2(y(i)), ...
+                     [abstol, reltol]);
+    endfor
+endfunction
+
+## method tests
+%!test
+%! f = @(x, y, z) x .* y .* z;
+%! assert (integral3 (f, 0, 1, 0, 1, 0, 1), 0.125, 1e-10);
+%! assert (integral3 (f, 0, 1, 0, 1, 0, 1, "method", "tiled"), 0.125, 1e-10);
+%! assert (integral3 (f, 0, 1, 0, 1, 0, 1, "method", "iterated"), 0.125, 1e-10);
+%! assert (integral3 (f, 0, 1, 0, 1, 0, 1, "method", "auto"), 0.125, 1e-10);
+
+## vectorized = false test
+%!test
+%! f = @(x, y, z) x * y * z;
+%! assert (integral3 (f, 0, 1, 0, 1, 0, 1, "vectorized", false), 0.125, 1e-10);
+
+## tolerance tests
+%!test
+%! f = @(x, y, z) 2 * x.^2 + 3 * y.^2 + 4 * z.^2;
+%! assert (integral3 (f, 0, 5, -5, 0, 0, 5, "AbsTol", 1e-9), 9375, 1e-9);
+%! assert (integral3 (f, 0, 5, -5, 0, 0, 5, "RelTol", 1e-6), 9375, -1e-6);
+%! assert (integral3 (f, 0, 5, -5, 0, 0, 5, "RelTol", 1e-6, "AbsTol", 1e-9),
+%!         9375, 1e-9);
+
+## non rectangular region
+## This test is too slow with "iterated" method
+%!assert (integral3 (@(x,y,z) 1 ./ (x + y + z), 0, 1, 0, @(x) 1 - x, 0, 
+%!        @(x, y) 1 - x - y, "method", "tiled"), 0.25, 1e-6)
+
+## Test input validation
+%!error integral3
+%!error integral3 (0, 1 ,2 ,3 ,4, 5, 6)
+%!error integral3 (@plus)
+%!error integral3 (@plus, 1)
+%!error integral3 (@plus, 1, 2)
+%!error integral3 (@plus, 1, 2, 3)
+%!error integral3 (@plus, 1, 2, 3, 4, 5, [6 7])
+%!error integral3 (@plus, 1, 2, 3, 4, 5, "test")
+%!error integral3 (@plus, 1, 2, 3, 4, 5, 6, "foo")
+%!error <unknown property 'foo'>  integral3 (@plus, 1, 2, 3, 4, 6, 6, "foo", "bar")
+%!error <property PROP must be a string> integral3 (@plus, 1, 2, 3, 4, 5, 6, 99, "bar")
+%!error <AbsTol value must be a numeric> integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", "foo")
+%!error <AbsTol value must be a .* scalar> integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", [1, 2])
+%!error <AbsTol value must be.* .= 0> integral3 (@plus, 1, 2, 3, 4, 5, 6, "AbsTol", -1)
+%!error <RelTol value must be a numeric> integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", "foo")
+%!error <RelTol value must be a .* scalar> integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", [1, 2])
+%!error <RelTol value must be.* .= 0> integral3 (@plus, 1, 2, 3, 4, 5, 6, "RelTol", -1)
+%!error <method 'bad' unrecognized> integral3 (@plus, 1, 2, 3, 4, 5, 6, "method", "bad")
diff -r d72ba9191352 -r 2c820ab8a0fb scripts/general/module.mk
--- a/scripts/general/module.mk	Wed Oct 04 14:13:49 2017 -0500
+++ b/scripts/general/module.mk	Sat Oct 14 18:30:51 2017 -0400
@@ -66,6 +66,7 @@
   %reldir%/postpad.m \
   %reldir%/prepad.m \
   %reldir%/publish.m \
+  %reldir%/quad2d.m \
   %reldir%/quadgk.m \
   %reldir%/quadl.m \
   %reldir%/quadv.m \
diff -r d72ba9191352 -r 2c820ab8a0fb scripts/general/quad2d.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/quad2d.m	Sat Oct 14 18:30:51 2017 -0400
@@ -0,0 +1,427 @@
+## Copyright (C) 2017 David Bateman
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {} {@var{q} =} quad2d (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb})
+## @deftypefnx {} {@var{q} =} quad2d (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{prop}, @var{val}, @dots{})
+## @deftypefnx {} {[@var{q}, @var{err}, @var{iter}] =} quad2d (@dots{})
+##
+## Numerically evaluate the two-dimensional integral of @var{f} using adaptive
+## quadrature over the two-dimensional domain defined by @var{xa}, @var{xb},
+## @var{ya}, @var{yb} using tiled integration. Additionally, @var{ya} and 
+## @var{yb} may be scalar functions of @var{x}, allowing for the integration 
+## over non-rectangular domains.
+##
+## @var{f} is a function handle, inline function, or string containing the name
+## of the function to evaluate.  The function @var{f} must be of the form
+## @math{z = f(x,y)} where @var{x} is a vector and @var{y} is a scalar.  It
+## should return a vector of the same length and orientation as @var{x}.
+##
+## Additional optional parameters can be specified using
+## @qcode{"@var{property}", @var{value}} pairs.  Valid properties are:
+##
+## @table @code
+## @item MaxFunEvals
+## The maximum number of function calls to the vectorized function @var{f}.
+## The default value is 5000. 
+##
+## @item FailurePlot
+## If @code{quad2d} fails to converge to the desired error tolerance before 
+## @var{MaxFunEvals} is reached, a plot of the areas that still need refinement 
+## is created. The default value is @var{false}.
+##
+## @item Singular
+## Use transforms to weaken singularities on the edge of the integration
+## doamin. Valid options are @var{true} or @var{false}. The default value is 
+## @var{true}.
+##
+## @item Vectorized
+## Option to disable vectorized integration, forcing octave to use only scalar 
+## inputs when calling the integrand. The default value is @var{false}.
+##
+## @item AbsTol
+## Define the absolute error tolerance for the quadrature.  The default
+## value is 1e-10 (1e-5 for single).
+##
+## @item RelTol
+## Define the relative error tolerance for the quadrature.  The default
+## value is 1e-6 (1e-4 for single).
+## @end table
+##
+## Adaptive quadrature is used to minimize the estimate of error until the
+## following is satisfied:
+## @tex
+## $$error \leq \max \left( AbsTol, RelTol\cdot\vert q\vert \right)$$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+##   @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|).
+## @end group
+## @end example
+##
+## @end ifnottex
+##
+## The optional output @var{err} is an approximate bound on the error in the 
+## integral @code{abs (@var{q} - @var{I})}, where @var{I} is the exact value 
+## of the integral, and @var{iter} is the number of vectorized function calls 
+## to the function @var{f} that were used.
+##
+## Known @sc{matlab} incompatibilities:
+##
+## @enumerate
+## @item
+## If tolerances are left unspecified, and any integration limits are
+## of type @code{single}, then Octave's integral functions automatically
+## reduce the default absolute and relative error tolerances as specified
+## above.  If tighter tolerances are desired they must be specified.
+## @sc{matlab} leaves the tighter tolerances appropriate for @code{double}
+## inputs in place regardless of the class of the integration limits.
+## @end enumerate
+##
+## Reference: @nospell{L.F. Shampine},
+## @cite{"@sc{matlab} program for quadrature in 2D"}, Applied Mathematics and 
+## Computation, pp. 266--274, Vol 1, 2008.
+##
+## @seealso{integral, quad2d, integral3, quad, quadgk, quadv, quadl, quadcc,
+##          trapz, dblquad, triplequad}
+## @end deftypefn
+
+function [q, err, iter] = quad2d (f, xa, xb, ya, yb, varargin)
+  
+  if ((nargin < 5) || (mod (nargin, 2) == 0))
+    print_usage ();
+  endif
+
+  ## Convert function given as a string to a function handle
+  if (ischar (f))
+    f = @(x) feval (f, x);
+  endif
+  
+  if (! is_function_handle (f))
+    print_usage ();
+  endif
+
+  if (! ((isscalar (xa)) && (isscalar (xb))))
+    print_usage ();
+  endif
+
+  ## Check for single or double limits to set appropriate default tolerance.
+  issingle = isa ([xa, xb], "single");
+  issingle = issingle && (! is_function_handle (ya)) && (isa (ya, "single"));
+  issingle = issingle && (! is_function_handle (yb)) && (isa (yb, "single"));
+
+  ## Set defaults, update with any specified parameters.
+  if issingle
+    abstol = 1e-5;
+    reltol = 1e-4;
+  else
+    abstol = 1e-10;
+    reltol = 1e-6;
+  endif
+  
+  vectorized = true;
+  singular = true;
+  idx = 1;
+  maxiter = 5000;
+  failureplot = false;
+  
+  while (idx < nargin - 5)
+    prop = varargin{idx++};
+    if (! ischar (prop))
+      error ("quad2d: property PROP must be a string");
+    endif
+
+    switch (tolower (prop))
+      case "abstol"
+        abstol = varargin{idx++};
+        if (! ((isnumeric (abstol)) && (isscalar (abstol)) && (abstol >= 0)))
+          error ("quad2d: AbsTol value must be a numeric scalar >= 0");
+        endif
+
+      case "reltol"
+        reltol = varargin{idx++};
+        if (! ((isnumeric (reltol)) && (isscalar (reltol)) && (reltol >= 0)))
+          error ("quad2d: RelTol value must be a numeric scalar >= 0");
+        endif
+     
+      case "vectorized"
+        vectorized = varargin{idx++};
+        if (! islogical (vectorized))
+          error ("quad2d: 'vectorized' must be a logical value");
+        endif
+        
+      case "singular"
+        singular = varargin{idx++};
+        if (! islogical (singular))
+          error ("quad2d: 'singular' must be a logical value");
+        endif
+        
+      case "maxfunevals"
+        maxiter = varargin{idx++};
+        if (! ((isnumeric (maxiter)) && (isscalar (maxiter)) && 
+            (floor (maxiter) == maxiter) && (maxiter >= 1)))
+          error ("quad2d: MaxFunEvals value must be an integer scalar >= 1");
+        endif
+        
+      case "failureplot"
+        failureplot = varargin{idx++};
+        if (! islogical (failureplot))
+          error ("quad2d: 'failureplot' must be a logical value");
+        endif
+      
+      otherwise
+        error ("quad2d: unknown property '%s'", prop);
+    endswitch
+  endwhile
+
+  if (! vectorized)
+    f = @(x, y) arrayfun (f, x, y);
+  endif
+  
+  # check upper and lower bounds of y
+  if (! is_function_handle (ya))
+    if (isscalar (ya))
+      ya = @(x) ya * ones(rows (x), columns (x));
+    else
+      error ("quad2d: 'ya' must be a constant or a (vectorized) function.");
+    endif
+  endif
+  if (! is_function_handle (yb))
+    if (isscalar (yb))
+      yb = @(x) yb * ones(rows (x), columns (x));
+    else
+      error ("quad2d: 'ya' must be a constant or a (vectorized) function.");
+    endif
+  endif
+
+    iter = 0;  
+    qaccept = 0;
+    qerraccept = 0;
+
+    if (singular)    
+      # Shampine suggests using the singularity weakening transform 
+      # suggested by Havie
+      #   \int_a^b f(x) dx = \int_0^pi f (g(t)) (dx / dt) dt
+      # where
+      #   g(t) = ((a - b) * cos(t) + (a + b)) / 2
+      #   dx = - (a - b) * sin(t) / 2 dt
+      # Now our integral is
+      #   \int_a^b \int_0^1 f(x,y) dydx
+      # as we already subsitute for "y", so
+      #   gx(tx) = ((a - b) * cos(tx) + (a + b)) / 2
+      #   gy(ty) = (1 - cos(ty)) / 2
+      #   dydx = (b - a) * sin(tx) * sin(ty) / 4 dtydtx
+      
+      xtrans = @(tx) ((xa - xb) .* cos (tx) + (xa + xb)) ./ 2;
+      ytrans = @(ty) (1 - cos (ty)) ./ 2;
+      ztrans = @(tx, ty) (xb - xa) .* sin(tx) .* sin(ty) ./ 4; 
+      area = pi ^ 2;
+
+      # Initialize tile list
+      tilelist(1) = struct("xa", 0, "xb", pi, "ya", 0, "yb", pi, "q", 0, ...
+                            "qerr", Inf);    
+    else
+      xtrans = @(tx) tx;
+      ytrans = @(ty) ty;
+      ztrans = @(tx, ty)  1;
+      area = (xb - xa);
+
+      # Initialize tile list
+      tilelist(1) = struct("xa", xa, "xb", xb, "ya", 0, "yb", 1, "q", 0, ...
+                            "qerr", Inf);
+    endif
+        
+    while ((length (tilelist) > 0) && (iter < maxiter))
+      # Get tile with the largest error
+      [~, indx] = max ([tilelist.qerr]);
+      tile = tilelist(indx);
+      tilelist(indx) = [];
+      
+      # Subdivide the tile into 4 subtiles
+      iter += 4;
+      tiles(4) = struct("xa", tile.xa, "xb", tile.xa + (tile.xb - tile.xa) / 2, 
+                        "ya", tile.ya, "yb", tile.ya + (tile.yb - tile.ya) / 2,
+                        "q", 0, "qerr", 0);
+      tiles(3) = struct("xa", tile.xa, "xb", tile.xa + (tile.xb - tile.xa) / 2, 
+                        "ya", tile.ya + (tile.yb - tile.ya) / 2, "yb", tile.yb,
+                        "q", 0, "qerr", 0);
+      tiles(2) = struct("xa", tile.xa + (tile.xb - tile.xa) / 2, "xb", tile.xb, 
+                        "ya", tile.ya, "yb", tile.ya + (tile.yb - tile.ya) / 2,
+                        "q", 0, "qerr", 0);
+      tiles(1) = struct("xa", tile.xa + (tile.xb - tile.xa) / 2, "xb", tile.xb, 
+                        "ya", tile.ya + (tile.yb - tile.ya) / 2, "yb", tile.yb,
+                        "q", 0, "qerr", 0);
+
+      # Perform the quadrature of 4 subtiles
+      for i = 1:4
+        [tiles(i).q, tiles(i).qerr] = tensorproduct (f, ya, yb, tiles(i), ...
+                                            xtrans, ytrans, ztrans, singular);
+      endfor
+ 
+      q = qaccept + sum ([[tilelist.q], tiles.q]);
+      err = qerraccept + sum ([[tilelist.qerr], tiles.qerr]);
+      tol = max (abstol, reltol .* abs (q));
+
+      # Shampine suggests taking a margin of a factor of 8 for 
+      # the local tolerance. That and the fact that we are subdividing
+      # into 4 tiles, means we divide by 32 at this point.
+      localtol = tol * ([tile.xb] - [tile.xa]) * ([tile.yb] - [tile.ya]) ...
+                / area / 32;
+      
+      # If global tolerance is met, return
+      if (err < tol)
+        break;
+      endif
+  
+      # Accept the tiles meeting the tolerance, and add the others back to
+      # the list of tiles to treat
+      indx = find ([tiles.qerr] < localtol);
+      qaccept += sum ([tiles(indx).q]);
+      qerraccept += sum ([tiles(indx).qerr]);
+      tiles(indx) = [];
+      tilelist = [tilelist, tiles];
+    endwhile
+    
+    if (iter >= maxiter)
+      if (err > max (abstol, reltol .* abs (q)))
+        warning ("quad2d : Maximum number of sub-tiles without convergence.");
+      else
+        warning(["quad2d : Maximum number of sub-tiles reached, accuracy ", ...
+                 "may be low."]);
+      endif
+      if (failureplot)
+        newplot;
+        title("QUAD2D : Areas needing refinement");
+        for tile = tilelist
+          xaa = xtrans(tile.xa);
+          xbb = xtrans(tile.xb);
+          y1 = ya(xaa) + ytrans(tile.ya) * (yb(xaa) - ya(xaa));
+          y2 = ya(xaa) + ytrans(tile.yb) * (yb(xaa) - ya(xaa));
+          y3 = ya(xbb) + ytrans(tile.yb) * (yb(xbb) - ya(xbb));
+          y4 = ya(xbb) + ytrans(tile.ya) * (yb(xbb) - ya(xbb));
+          patch ([xaa, xaa, xbb, xbb, xaa],[y1, y2, y3, y4, y1], 'b');
+        endfor
+      endif
+    endif
+endfunction
+
+function [q, qerr] = tensorproduct (f, ya, yb, tile, xtrans, ytrans, ...
+                                     ztrans, singular)
+  # The Shampine TwoD paper proposes using a G3,K7 rule in a tensor product
+  # I couldn't find a tabulated abscissas and weights of a G3,K7 rule publically
+  # available, so I'm using a G7,K15 rule from Octave's implementation of quadgk 
+  
+  persistent abscissa = [-0.9914553711208126e+00, -0.9491079123427585e+00, ...
+                         -0.8648644233597691e+00, -0.7415311855993944e+00, ...
+                         -0.5860872354676911e+00, -0.4058451513773972e+00, ...
+                         -0.2077849550078985e+00,  0.0000000000000000e+00, ...
+                          0.2077849550078985e+00,  0.4058451513773972e+00, ...
+                          0.5860872354676911e+00,  0.7415311855993944e+00, ...
+                          0.8648644233597691e+00,  0.9491079123427585e+00, ...
+                          0.9914553711208126e+00];
+
+  persistent weights15 = ...
+            [0.2293532201052922e-01,  0.6309209262997855e-01, ...
+             0.1047900103222502e+00,  0.1406532597155259e+00, ...
+             0.1690047266392679e+00,  0.1903505780647854e+00, ...
+             0.2044329400752989e+00,  0.2094821410847278e+00, ...
+             0.2044329400752989e+00,  0.1903505780647854e+00, ...
+             0.1690047266392679e+00,  0.1406532597155259e+00, ...
+             0.1047900103222502e+00,  0.6309209262997855e-01, ...
+             0.2293532201052922e-01];
+
+  persistent weights7  = [0.0, ...
+             0.1294849661688697e+00, 0.0, ...
+             0.2797053914892767e+00, 0.0, ...
+             0.3818300505051889e+00, 0.0, ...
+             0.4179591836734694e+00, 0.0, ...
+             0.3818300505051889e+00, 0.0, ...
+             0.2797053914892767e+00, 0.0, ...
+             0.1294849661688697e+00, 0.0];
+  
+  xaa = tile.xa;
+  xbb = tile.xb;
+  yaa = tile.ya;
+  ybb = tile.yb;
+  
+  tx = ((xbb - xaa) * abscissa + xaa + xbb) / 2;
+  x = xtrans(tx);
+  ty = (abscissa' * (ybb - yaa)  + yaa + ybb) / 2;
+  y = ones (15, 1) * ya(x) + ytrans(ty) * (yb(x) - ya(x));
+  
+  xhalfwidth = (xbb - xaa ) / 2;
+  yhalfwidth = ones (15, 1) * (yb(x) - ya(x)) .* (ybb - yaa) ./ 2;
+  
+  x = ones (15, 1) * x;  
+  tx = ones (15,1) * tx;
+  ty = ty * ones (1, 15);
+  
+  z = yhalfwidth .* f (x, y) .* ztrans(tx, ty) .* xhalfwidth;
+  q = weights15 * (weights15 * z)';
+  qerr = abs(weights7 * (weights7 * z)' - q);
+endfunction
+
+## Test input validation
+%!error quad2d
+%!error quad2d (0, 1 ,2 ,3 ,4)
+%!error quad2d (@plus)
+%!error quad2d (@plus, 1)
+%!error quad2d (@plus, 1, 2)
+%!error quad2d (@plus, 1, 2, 3)
+%!error quad2d (@plus, 1, 2, 3, [4 5])
+%!error quad2d (@plus, 1, 2, 3, "test")
+%!error quad2d (@plus, 1, 2, 3, 4, "foo")
+%!error <unknown property 'foo'>  quad2d (@plus, 1, 2, 3, 4, "foo", "bar")
+%!error <property PROP must be a string> quad2d (@plus, 1, 2, 3, 4, 99, "bar")
+%!error <AbsTol value must be a numeric> quad2d (@plus, 1, 2, 3, 4, "AbsTol", "foo")
+%!error <AbsTol value must be a .* scalar> quad2d (@plus, 1, 2, 3, 4, "AbsTol", [1, 2])
+%!error <AbsTol value must be.* .= 0> quad2d (@plus, 1, 2, 3, 4, "AbsTol", -1)
+%!error <RelTol value must be a numeric> quad2d (@plus, 1, 2, 3, 4, "RelTol", "foo")
+%!error <RelTol value must be a .* scalar> quad2d (@plus, 1, 2, 3, 4, "RelTol", [1, 2])
+%!error <RelTol value must be.* .= 0> quad2d (@plus, 1, 2, 3, 4, "RelTol", -1)
+%!warning <Maximum number of> quad2d (@plus, 1, 2, 3, 4, "MaxFunEvals", 1);
+
+%!test
+%! f = @(x, y) x .* y;
+%! assert (quad2d (f, 0, 1, 0, 1), 0.25, 1e-10);
+
+%!test
+%! f = @(x, y) 9 * x.^2 + 15 * y.^2;
+%! assert (quad2d (f, 0, 5, -5, 0, "AbsTol", 1e-9), 5000, 1e-9);
+%! assert (quad2d (f, 0, 5, -5, 0, "RelTol", 1e-6), 5000, -1e-6);
+%! assert (quad2d (f, 0, 5, -5, 0, "RelTol", 1e-6, "AbsTol", 1e-9),
+%!         5000, 1e-9);
+
+## tests from dblquad
+%!assert (quad2d (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-7),
+%!        2*log (2), 1e-6)
+%!assert (quad2d (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "RelTol", 1e-6),
+%!        2*log (2), -1e-6)
+%!assert (quad2d (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-8,
+%!        "RelTol", 1e-6), 2*log (2), -1e-6)
+%!assert (quad2d (@(x, y) exp (-x.^2 - y.^2) , -1, 1, -1, 1),
+%!        pi * erf (1).^2, 1e-10)
+
+%!assert (quad2d (@plus, 1, 2, 3, 4), 5, 1e-10)
+
+%!assert (quad2d (@(x,y) 1 ./ (x + y), 0, 1, 0, @(x) 1 - x), 1, -1e-6)
+
+%!assert (quad2d (@(x,y) 1 ./ (x + y), 0, 1, 0, @(x) 1 - x, "Singular", true),
+%!        1, -1e-6)
\ No newline at end of file
diff -r d72ba9191352 -r 2c820ab8a0fb scripts/help/__unimplemented__.m
--- a/scripts/help/__unimplemented__.m	Wed Oct 04 14:13:49 2017 -0500
+++ b/scripts/help/__unimplemented__.m	Sat Oct 14 18:30:51 2017 -0400
@@ -66,7 +66,7 @@
 
     case "matlabrc"
       txt = ["matlabrc is not implemented.  ", ...
-             'Octave uses the file ".octaverc" instead.'];
+             "Octave uses the file '.octaverc' instead."];
 
     case {"ode113", "ode15i", "ode15s", "ode23s", "ode23t", "ode23tb"}
       txt = ["Octave provides lsode and ode45 for solving differential equations. ", ...
@@ -78,9 +78,6 @@
       txt = ["startup is not implemented.  ", ...
              'Octave uses the file ".octaverc" instead.'];
 
-    case "quad2d"
-      txt = ["quad2d is not implemented.  Consider using dblquad."];
-
     case {"xlsread", "xlsfinfo", "xlswrite", "wk1read", "wk1finfo", "wk1write"}
       txt = ["Functions for spreadsheet style I/O ", ...
              "(.xls .xlsx .sxc .ods .dbf .wk1 etc.) " , ...
@@ -843,11 +840,10 @@
   };
 endfunction
 
-
 %!test
 %! str = __unimplemented__ ("no_name_function");
 %! assert (isempty (str));
-%! str = __unimplemented__ ("quad2d");
-%! assert (str(1:51), "quad2d is not implemented.  Consider using dblquad.");
+%! str = __unimplemented__ ("matlabrc");
+%! assert (str(1:71), "matlabrc is not implemented.  Octave uses the file '.octaverc' instead.");
 %! str = __unimplemented__ ("MException");
 %! assert (str(1:58), "the 'MException' function is not yet implemented in Octave");
