--- a/NEWS	Mon Sep 25 14:59:19 2017 -0700
+++ b/NEWS	Sat Sep 23 18:11:20 2017 -0400
@@ -67,6 +67,8 @@
       hgtransform
       humps
       integral
+      integral2
+      integral3
       openvar
       repelem
       rticks

--- a/doc/interpreter/quad.txi	Mon Sep 25 14:59:19 2017 -0700
+++ b/doc/interpreter/quad.txi	Sat Sep 23 18:11:20 2017 -0400
@@ -319,8 +319,8 @@
 @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:
+@code{triplequad} functions for integrals over two and three variables.  For
+example:
 
 @example
 @group
@@ -333,6 +333,13 @@
 
 @DOCSTRING(triplequad)
 
+@code{integral2} and @code{integral3} are also provided as partial
+compatibility wrappers to @code{dblquad} and @code{triplequad}:
+
+@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

--- a/scripts/general/dblquad.m	Mon Sep 25 14:59:19 2017 -0700
+++ b/scripts/general/dblquad.m	Sat Sep 23 18:11:20 2017 -0400
@@ -42,7 +42,8 @@
 ## Additional arguments, are passed directly to @var{f}.  To use the default
 ## value for @var{tol} or @var{quadf} one may pass @qcode{':'} or an empty
 ## matrix ([]).
-## @seealso{triplequad, quad, quadv, quadl, quadgk, quadcc, trapz}
+## @seealso{integral2, integral3, triplequad, quad, quadv, quadl, quadgk,
+##          quadcc, trapz}
 ## @end deftypefn
 
 function q = dblquad (f, xa, xb, ya, yb, tol = 1e-6, quadf = @quadcc, varargin)

--- a/scripts/general/integral.m	Mon Sep 25 14:59:19 2017 -0700
+++ b/scripts/general/integral.m	Sat Sep 23 18:11:20 2017 -0400
@@ -24,8 +24,8 @@
 ## adaptive quadrature.
 ##
 ## @code{integral} is a wrapper for @code{quadgk} (for scalar integrands) and
-## @code{quadv} (for array-valued integrands) intended to provide Matlab
-## compatibility. More control of the numerical integration may be achievable
+## @code{quadv} (for array-valued integrands) 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
@@ -36,20 +36,12 @@
 ## or both limits may be infinite or contain weak end singularities.  If either
 ## or both limits are complex, @code{integral} will perform a straight line
 ## path integral.  Alternatively, a complex domain path can be specified using
-## the "Waypoints" option (see below).
+## the @qcode{"Waypoints"} option (see below).
 ##
 ## 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
-## absolute tolerance is 1e-10 (1e-5 for single).
-##
-## @item RelTol
-## Define the relative error tolerance for the quadrature.  The default
-## relative tolerance is 1e-6 (1e-4 for single).
-##
 ## @item Waypoints
 ## Specifies points to be used in defining subintervals of the quadrature
 ## algorithm, or if @var{a}, @var{b}, or @var{waypoints} are complex then
@@ -61,14 +53,50 @@
 ## @var{arrayvalued} is specified as true.  This option will cause
 ## @code{integral} to perform the integration over the entire array and return
 ## @var{q} with the same dimensions as returned by @var{f}.
+##
+## @item AbsTol
+## Define the absolute error tolerance for the quadrature.  The default
+## absolute tolerance is 1e-10 (1e-5 for single).
+##
+## @item RelTol
+## Define the relative error tolerance for the quadrature.  The default
+## relative tolerance is 1e-6 (1e-4 for single).
 ## @end table
 ##
-## Implementation Note: As a consequence of using @code{quadgk} and
-## @code{quadv}, certain option combinations are currently unsupported.
-## @qcode{"ArrayValued"} cannot be combined with @qcode{"RelTol"} or
-## @qcode{"Waypoints"}.  This is a known incompatibility with Matlab.
+## 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:
 ##
-## @seealso{quad, quadgk, quadv, quadl, quadcc, trapz, dblquad, triplequad}
+## @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
+## As a consequence of using @code{quadgk} and @code{quadv}, certain option
+## combinations are not supported.  Currently, @qcode{"ArrayValued"} cannot be
+## combined with @qcode{"RelTol"} or @qcode{"Waypoints"}.
+## @end enumerate
+##
+## @seealso{integral2, integral3, quad, quadgk, quadv, quadl, quadcc, trapz,
+##          dblquad, triplequad}
 ## @end deftypefn
 
 function q = integral (f, a, b, varargin)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/integral2.m	Sat Sep 23 18:11:20 2017 -0400
@@ -0,0 +1,221 @@
+## 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");

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/integral3.m	Sat Sep 23 18:11:20 2017 -0400
@@ -0,0 +1,217 @@
+## 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");

--- a/scripts/general/module.mk	Mon Sep 25 14:59:19 2017 -0700
+++ b/scripts/general/module.mk	Sat Sep 23 18:11:20 2017 -0400
@@ -42,6 +42,8 @@
   %reldir%/inputParser.m \
   %reldir%/int2str.m \
   %reldir%/integral.m \
+  %reldir%/integral2.m \
+  %reldir%/integral3.m \
   %reldir%/interp1.m \
   %reldir%/interp2.m \
   %reldir%/interp3.m \

--- a/scripts/general/quadgk.m	Mon Sep 25 14:59:19 2017 -0700
+++ b/scripts/general/quadgk.m	Sat Sep 23 18:11:20 2017 -0400
@@ -45,7 +45,7 @@
 ## same.
 ##
 ## The optional argument @var{tol} defines the absolute tolerance used to stop
-## the integration procedure.  The default value is 1e-10.
+## the integration procedure.  The default value is 1e-10 (1e-5 for single).
 ##
 ## The algorithm used by @code{quadgk} involves subdividing the integration
 ## interval and evaluating each subinterval.  If @var{trace} is true then after
@@ -120,7 +120,8 @@
 ## Computational and Applied Mathematics, pp. 131--140, Vol 211, Issue 2,
 ## Feb 2008.
 ##
-## @seealso{quad, quadv, quadl, quadcc, trapz, dblquad, triplequad}
+## @seealso{quad, quadv, quadl, quadcc, trapz, dblquad, triplequad, integral, 
+##           integral2, integral3}
 ## @end deftypefn
 
 function [q, err] = quadgk (f, a, b, varargin)

--- a/scripts/general/quadl.m	Mon Sep 25 14:59:19 2017 -0700
+++ b/scripts/general/quadl.m	Sat Sep 23 18:11:20 2017 -0400
@@ -54,7 +54,8 @@
 ## Reference: @nospell{W. Gander and W. Gautschi}, @cite{Adaptive Quadrature -
 ## Revisited}, BIT Vol. 40, No. 1, March 2000, pp. 84--101.
 ## @url{http://www.inf.ethz.ch/personal/gander/}
-## @seealso{quad, quadv, quadgk, quadcc, trapz, dblquad, triplequad}
+## @seealso{quad, quadv, quadgk, quadcc, trapz, dblquad, triplequad, integral, 
+##          integral2, integral3}
 ## @end deftypefn
 
 ## Original Author: Walter Gautschi

--- a/scripts/general/quadv.m	Mon Sep 25 14:59:19 2017 -0700
+++ b/scripts/general/quadv.m	Sat Sep 23 18:11:20 2017 -0400
@@ -57,7 +57,8 @@
 ## Note: @code{quadv} is written in Octave's scripting language and can be
 ## used recursively in @code{dblquad} and @code{triplequad}, unlike the
 ## @code{quad} function.
-## @seealso{quad, quadl, quadgk, quadcc, trapz, dblquad, triplequad}
+## @seealso{quad, quadl, quadgk, quadcc, trapz, dblquad, triplequad, integral, 
+##          integral2, integral3}
 ## @end deftypefn
 
 function [q, nfun] = quadv (f, a, b, tol = [], trace = false, varargin)

--- a/scripts/general/triplequad.m	Mon Sep 25 14:59:19 2017 -0700
+++ b/scripts/general/triplequad.m	Sat Sep 23 18:11:20 2017 -0400
@@ -43,7 +43,8 @@
 ## Additional arguments, are passed directly to @var{f}.  To use the default
 ## value for @var{tol} or @var{quadf} one may pass @qcode{':'} or an empty
 ## matrix ([]).
-## @seealso{dblquad, quad, quadv, quadl, quadgk, quadcc, trapz}
+## @seealso{integral3, integral2, dblquad, quad, quadv, quadl, quadgk, quadcc,
+##          trapz}
 ## @end deftypefn
 
 function q = triplequad (f, xa, xb, ya, yb, za, zb, tol = 1e-6, quadf = @quadcc, varargin)

--- a/scripts/help/__unimplemented__.m	Mon Sep 25 14:59:19 2017 -0700
+++ b/scripts/help/__unimplemented__.m	Sat Sep 23 18:11:20 2017 -0400
@@ -59,12 +59,6 @@
       txt = ["griddedInterpolant is not implemented.  ", ...
              "Consider using griddata."];
 
-    case "integral2"
-      txt = ["integral2 is not implemented.  Consider using dblquad."];
-
-    case "integral3"
-      txt = ["integral3 is not implemented.  Consider using triplequad"];
-
     case "linprog"
       txt = ["Octave does not currently provide linprog.  ", ...
              "Linear programming problems may be solved using @code{glpk}.  ", ...