  The algorithm AMGKQ for adaptive multivariate Gauss-Kronrod quadrature over
hyper-rectangular regions of arbitrary dimensionality is proposed and
implemented in Octave/MATLAB. It can approximate numerically any number of
integrals over a common domain simultaneously. Improper integrals are addressed
through singularity weakening coordinate transformations. Internal
singularities are addressed through the use of breakpoints. Its accuracy
performance is investigated thoroughly, and its running time is compared to
other commonly available routines in two and three dimensions. Its running time
can be several orders of magnitude faster than recursively called quadrature
routines. Its performance is limited only by the memory structure of its
operating environment. Included with the software are numerous examples of its
invocation.
