--- a/scripts/general/integral2.m	Tue Jan 30 20:02:05 2024 -0500
+++ b/scripts/general/integral2.m	Wed Jan 31 00:21:31 2024 -0500
@@ -34,10 +34,17 @@
 ## @var{ya} and @var{yb} may be scalar functions of @var{x}, allowing for
 ## 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}.
+## @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)}, and all operations must be vectorized such that @var{x}
+## and @var{y} accept array inputs and return array outputs of the same size.
+## (It can be assumed that @var{x} and @var{y} will either be same-size arrays
+## or one will be a scalar.)  The underlying integrators will input arrays of
+## integration points into @var{f} and/or use internal vector expansions to
+## speed computation that can produce unpredictable results if @var{f} is not
+## restricted to elementwise operations.  For integrands where this is
+## unavoidable, the @qcode("Vectorized") option described below may produce
+## more reliable results.
 ##
 ## Additional optional parameters can be specified using
 ## @qcode{"@var{property}", @var{value}} pairs.  Valid properties are:
@@ -60,9 +67,10 @@
 ## @item Vectorized
 ## Enable or disable vectorized integration.  A value of @code{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 to be used.  The default value is
-## @code{true}.
+## integrands @math{f(x,y)} that have not been vectorized or only accept
+## scalar values of @var{x} or @var{y}.  The default value is @code{true}.
+## Note that this is achieved by wrapping @math{f(x,y)} with the function
+## @code{arrayfun}, which may significantly decrease computation speed.
 ## @end table
 ##
 ## Adaptive quadrature is used to minimize the estimate of error until the
@@ -95,7 +103,7 @@
 ## @end example
 ##
 ## The result is a volume, which for this constant-value integrand, is just
-## @code{@var{Length} * @var{Width} * @var{Height}}.
+## @w{@code{@var{Length} * @var{Width} * @var{Height}}}.
 ##
 ## Example 2 : integrate a triangular region in x-y plane
 ##
@@ -108,9 +116,34 @@
 ## @end group
 ## @end example
 ##
-## The result is a volume, which for this constant-value integrand, is the
-## Triangle Area x Height or
-## @code{1/2 * @var{Base} * @var{Width} * @var{Height}}.
+## The result is a volume, which for this constant-value integrand
+## @w{@math{@var{f} = 2}}, is the Triangle Area x Height or
+## @w{@code{1/2 * @var{Base} * @var{Width} * @var{Height}}}.
+##
+## Example 3 : integrate a non-vectorized function over a square region
+##
+## @example
+## @group
+## @var{f} = @@(@var{x},@var{y}) sinc (@var{x}) * sinc (@var{y}));
+## @var{q} = integral2 (@var{f}, -1, 1, -1, 1)
+##   @result{} @var{q} =  12.328  (incorrect)
+## @var{q} = integral2 (@var{f}, -1, 1, -1, 1, "Vectorized", false)
+##   @result{} @var{q} =  1.390 (correct)
+## @var{f} = @@(@var{x},@var{y}) sinc (@var{x}) .* sinc (@var{y});
+## @var{q} = integral2 (@var{f}, -1, 1, -1, 1)
+##   @result{} @var{q} =  1.390  (correct)
+## @end group
+## @end example
+##
+## The first result is incorrect as the non-elementwise operator between the
+## sinc functions in @var{f} create unintended matrix multiplications between
+## the internal integration arrays used by @code{integral2}.  In the second
+## result, setting @qcode{"Vectorized"} to false forces @code{integral2} to
+## perform scalar internal operations to compute the integral, resulting in
+## the correct numerical result at the cost of about a 20x increase in
+## computation time. In the third result, vectorizing the integrand @var{f}
+## using the elementwise multiplication operator gets the correct result
+## without increasing computation time.
 ##
 ## Programming Notes: If there are singularities within the integration region
 ## it is best to split the integral and place the singularities on the

--- a/scripts/general/integral3.m	Tue Jan 30 20:02:05 2024 -0500
+++ b/scripts/general/integral3.m	Wed Jan 31 00:21:31 2024 -0500
@@ -35,10 +35,17 @@
 ## functions of @var{x} and @var{y}, allowing for 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}.
+## @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,z)}, and all operations must be vectorized such that
+## @var{x}, @var{y}, and @var{z} accept array inputs and return array outputs
+## of the same size.  (It can be assumed that @var{x}, @var{y}, and @var{z}
+## will either be same-size arrays or scalars.)  The underlying integrators
+## will input arrays of integration points into @var{f} and/or use internal
+## vector expansions to speed computation that can produce unpredictable
+## results if @var{f} is not restricted to elementwise operations.  For
+## integrands where this is unavoidable, the @qcode("Vectorized") option
+## described below may produce more reliable results.
 ##
 ## Additional optional parameters can be specified using
 ## @qcode{"@var{property}", @var{value}} pairs.  Valid properties are:
@@ -61,9 +68,11 @@
 ## @item Vectorized
 ## Enable or disable vectorized integration.  A value of @code{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 to be used.  The default value is
-## @code{true}.
+## integrands @math{f(x,y,z)} that have not been vectorized or only accept
+## scalar values of @var{x}, @var{y}, or @var{z}.  The default value is
+## @code{true}.  Note that this is achieved by wrapping @math{f(x,y,z)} with
+## the function @code{arrayfun}, which may significantly decrease computation
+## speed.
 ## @end table
 ##
 ## Adaptive quadrature is used to minimize the estimate of error until the
@@ -113,6 +122,31 @@
 ## For this constant-value integrand, the result is a volume which is 1/8th
 ## of a unit sphere or @code{1/8 * 4/3 * pi}.
 ##
+## Example 3 : integrate a non-vectorized function over a cubic volume
+##
+## @example
+## @group
+## @var{f} = @@(@var{x},@var{y}) sinc (@var{x}) * sinc (@var{y}), * sinc (@var{z});
+## @var{q} = integral3 (@var{f}, -1, 1, -1, 1, -1, 1)
+##   @result{} @var{q} =  14.535  (incorrect)
+## @var{q} = integral3 (@var{f}, -1, 1, -1, 1, -1, 1, "Vectorized", false)
+##   @result{} @var{q} =  1.6388 (correct)
+## @var{f} = @@(@var{x},@var{y},@var{z}) sinc (@var{x}) .* sinc (@var{y}), .* sinc (@var{z});
+## @var{q} = integral3 (@var{f}, -1, 1, -1, 1, -1, 1)
+##   @result{} @var{q} =  1.6388  (correct)
+## @end group
+## @end example
+##
+## The first result is incorrect as the non-elementwise operator between the
+## sinc functions in @var{f} create unintended matrix multiplications between
+## the internal integration arrays used by @code{integral3}.  In the second
+## result, setting @qcode{"Vectorized"} to false forces @code{integral3} to
+## perform scalar internal operations to compute the integral, resulting in
+## the correct numerical result at the cost of about a 30x increase in
+## computation time. In the third result, vectorizing the integrand @var{f}
+## using the elementwise multiplication operator gets the correct result
+## without increasing computation time.
+##
 ## Programming Notes: If there are singularities within the integration region
 ## it is best to split the integral and place the singularities on the
 ## boundary.

--- a/scripts/general/quad2d.m	Tue Jan 30 20:02:05 2024 -0500
+++ b/scripts/general/quad2d.m	Wed Jan 31 00:21:31 2024 -0500
@@ -34,10 +34,17 @@
 ## @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}.
+## @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)}, and all operations must be vectorized such that @var{x}
+## and @var{y} accept array inputs and return array outputs of the same size.
+## (It can be assumed that @var{x} and @var{y} will either be same-size arrays
+## or one will be a scalar.)  The underlying integrators will input arrays of
+## integration points into @var{f} and/or use internal vector expansions to
+## speed computation that can produce unpredictable results if @var{f} is not
+## restricted to elementwise operations.  For integrands where this is
+## unavoidable, the @qcode("Vectorized") option described below may produce
+## more reliable results.
 ##
 ## Additional optional parameters can be specified using
 ## @qcode{"@var{property}", @var{value}} pairs.  Valid properties are:
@@ -60,8 +67,12 @@
 ## integration domain.  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{true}.
+## Enable or disable vectorized integration.  A value of @code{false} forces
+## Octave to use only scalar inputs when calling the integrand, which enables
+## integrands @math{f(x,y)} that have not been vectorized or only accept
+## scalar values of @var{x} or @var{y}.  The default value is @code{true}.
+## Note that this is achieved by wrapping @math{f(x,y)} with the function
+## @code{arrayfun}, which may significantly decrease computation speed.
 ##
 ## @item FailurePlot
 ## If @code{quad2d} fails to converge to the desired error tolerance before
@@ -113,9 +124,34 @@
 ## @end group
 ## @end example
 ##
-## The result is a volume, which for this constant-value integrand, is the
-## Triangle Area x Height or
-## @code{1/2 * @var{Base} * @var{Width} * @var{Height}}.
+## The result is a volume, which for this constant-value integrand
+## @w{@math{@var{f} = 2}}, is the Triangle Area x Height or
+## @w{@code{1/2 * @var{Base} * @var{Width} * @var{Height}}}.
+##
+## Example 3 : integrate a non-vectorized function over a square region
+##
+## @example
+## @group
+## @var{f} = @@(@var{x},@var{y}) sinc (@var{x}) * sinc (@var{y}));
+## @var{q} = quad2d (@var{f}, -1, 1, -1, 1)
+##   @result{} @var{q} =  12.328  (incorrect)
+## @var{q} = quad2d (@var{f}, -1, 1, -1, 1, "Vectorized", false)
+##   @result{} @var{q} =  1.390 (correct)
+## @var{f} = @@(@var{x},@var{y}) sinc (@var{x}) .* sinc (@var{y});
+## @var{q} = quad2d (@var{f}, -1, 1, -1, 1)
+##   @result{} @var{q} =  1.390  (correct)
+## @end group
+## @end example
+##
+## The first result is incorrect as the non-elementwise operator between the
+## sinc functions in @var{f} create unintended matrix multiplications between
+## the internal integration arrays used by @code{quad2d}.  In the second
+## result, setting @qcode{"Vectorized"} to false forces @code{quad2d} to
+## perform scalar internal operations to compute the integral, resulting in
+## the correct numerical result at the cost of about a 20x increase in
+## computation time. In the third result, vectorizing the integrand @var{f}
+## using the elementwise multiplication operator gets the correct result
+## without increasing computation time.
 ##
 ## Programming Notes: If there are singularities within the integration region
 ## it is best to split the integral and place the singularities on the