Using Legendre polynomials for the Gaussian quadrature is not very efficient. There are several reasons for this:

• You can easily end up in situations where the integrand diverges
• The limits \( \pm \infty \) have to be approximated with a finite number

It is very useful here to change to spherical coordinates $$ d{\bf r}_1d{\bf r}_2 = r_1^2dr_1 r_2^2dr_2 dcos(\theta_1)dcos(\theta_2)d\phi_1d\phi_2, $$ and $$ \frac{1}{r_{12}}= \frac{1}{\sqrt{r_1^2+r_2^2-2r_1r_2cos(\beta)}} $$ with $$ \cos(\beta) = \cos(\theta_1)\cos(\theta_2)+\sin(\theta_1)\sin(\theta_2)\cos(\phi_1-\phi_2)) $$