--- a/doc/interpreter/quad.txi
+++ b/doc/interpreter/quad.txi
@@ -50,6 +50,12 @@
 @item quadl
 Numerical integration using an adaptive Lobatto rule.
 
+@item quadgk
+Numerical integration using an adaptive Guass-Konrod rule.
+
+@item quadv
+Numerical integration using an adaptive vectorized Simpson's rule.
+
 @item trapz
 Numerical integration using the trapezoidal method.
 @end table
@@ -118,6 +124,10 @@
 
 @DOCSTRING(quadl)
 
+@DOCSTRING(quadgk)
+
+@DOCSTRING(quadv)
+
 @DOCSTRING(trapz)
 
 @DOCSTRING(cumtrapz)
@@ -173,10 +183,10 @@
 @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.  It is however possible to compute
-the integral of a function of multiple variables using the functions
-for one-dimensional integrals.
+Octave does not have built-in functions for computing the integral of
+functions of multiple variables directly.  It is however possible to
+compute the integral of a function of multiple variables using the
+functions for one-dimensional integrals.
 
 To illustrate how the integration can be performed, we will integrate
 the function
@@ -215,6 +225,19 @@
       @result{} 0.30022
 @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
+
+@example
+I =  dblquad (@(x, y) sin(pi.*x.*y).*sqrt(x.*y), 0, 1, 0, 1)
+      @result{} 0.30022
+@end example
+
+@DOCSTRING(dblquad)
+
+@DOCSTRING(triplequad)
+
 The above mentioned approach works but is fairly slow, and that problem
 increases exponentially with the dimensionality the problem.  Another
 possible solution is to use Orthogonal Collocation as described in the
@@ -250,4 +273,3 @@
 
 
 
-