Abstract: The algorithm AMGKQ for adaptive multivariate Gauss-Kronrod quadrature over hyperrectangular 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.