An optimization of cubature (OOC) method is developed to refine computationally intensive numerical integrations by optimally extending an existing cubature. Typically, to increase the accuracy of numerical integrations using methods which cannot be nested, such as Gaussian quadrature, the integrand is evaluated over a larger, disjoint set of abscissas, without using the previous integrand evaluations. The developed OOC method adds any number of abscissas to the existing quadrature and reevaluates all associated weights. To optimize the abscissas and weights, a global optimizationtechnique is used to minimize the sum of squared numerical integration error for a set of training functions, calculated as a function of the optimized weights and abscissas. The training functions must have a known integral over a hypercube domain. The abscissas are constrained to the domain and weights are constrained to be positive. The optimized weights and abscissas are then used to numerically integrate a function of over the domain. Additionally, a method for the optimization of weights, a variant of the OOC method in which only cubature weights are optimized, is presented and compared to traditional methods such as Gauss quadrature and Bayesian quadrature. The OOC method performs well on multivariate integrands, however, the optimization required by the OOC method may add significant computational expense. Therefore, the OOC method is determined to be best suited for computationally expensive, multivariate integrands for which the number of integrand evaluations is severely limited.