A globally adaptive algorithm for numerical cubature of a vector of functions over a collection of n -dimensional simplices is described. The algorithm is based on a subdivision strategy that chooses for subdivision at each stage the subregion (of the input simplices) with the largest estimated error. This subregion is divided into two, three or four equal volume subregions by cutting selected edges. These edges are selected using information about the smoothness of the integrands in the edge directions. The algorithm allows a choice from several embedded cubature rule sequences for approximate integration and error estimation. A Fortran 95 implementation as a part of CUBPACK is also discussed. Testing of the algorithm is described.