The authors address the problem of solving the system of nonlinear equations for the knots and weights of a cubature scheme for approximating an n-dimensional integral. The cubature formula is taken to be exact for all polynomials up to a certain degree, so the number of resulting equations is equal to the number of unknowns specifying the knots and weights. A continuation method is used to solve the resulting system of nonlinear equations in which the initial solution is obtained by causing the cubature formula to be exact for a class of piecewise affine functions on a triangulation of the domain of integration. The authors illustrate the technique for computation of integrals over a square by use of invariant (with respect to symmetries of the square) quadrature formulae.