In this article, we investigate an integral result presented by Herrick and Stillinger in 1975 [Phys. Rev. A 11, 42 (1975)] numerically and analytically. Based on their result we implement a custom integration routine. We demonstrate that custom made integration routines can be several times faster than the native routines. This may be beneficial in cases where many integrals need to be carried out. We compared the results of custom integration routine against those of the native integration routine. For cases, where the native routine could not perform the integration, we compared the custom integration routine against a Monte Carlo estimate. In all cases, we found excellent agreement. As an application of the implemented integration routine, we compute the $N$-body matrix elements of a Hamiltonian with two-body potentials.