I'm trying to calculate the total energy of a simple two charge system through the integral for electrostatic energy of a system given in Griffiths' book:
$$U = \frac{\epsilon_0}{2}\int_V E^2 dV .$$
Where the volume is integrated across all space so the boundary term not shown here decays to zero. I think that this should yield the same answer as the standard formula given for point charges:
$$U = \frac{1}{4\pi\varepsilon_0}\frac{Q_1Q_2}{R}.$$
But I'm having trouble evaluating the integral itself. I placed $Q_1$ on the origin of the coordinate axes and $Q_2$ on the $z$-axis a distance $R$ away from the first charge, and expanded the $E^2$ term:
$$E = E_1 + E_2 $$ so $$E^2 = E_1^2 + 2E_1 \centerdot E_2 + E_2^2.$$
I found that the integral of the self terms diverges when evaluated, and, after reading through Griffiths, decided to discard the self-energy terms and only retain the energy due to the exchange term.
Letting $r = \sqrt{x^2+y^2+z^2}$ and $r'= \sqrt{x^2+y^2+(z-R)^2}$, I found the integral of the interaction term to be:
$$E_1 = \frac{1}{4\pi\varepsilon_0}\frac{Q_1}{r^3}\vec{r}\quad\text{and}\quad E_2 \frac{1}{4\pi\varepsilon_0}\frac{Q_2}{r'^3}\vec{r'}$$
$$U = \epsilon_0\int_V E_1\centerdot E_2 \space dV = \frac{Q_1 Q_2}{16\pi^2\varepsilon_0}\int_V \frac{x^2 + y^2 + z^2-zR}{(x^2 + y^2 + z^2)^{\frac{3}{2}} \space (x^2+y^2+(z-R)^2)^{\frac{3}{2}}}\space dV.$$
Converting to spherical coordinates, with $r=\sqrt{x^2+y^2+z^2}$, $\theta $ the angle from the z-axis and $\varphi$ the azimutal angle, where I have evaluated the azimuthal integral:
$$U = \frac{Q_1 Q_2}{8\pi\varepsilon_0}\int_0^\infty \int_0^{2\pi} \frac{r - R\cos(\theta)}{(r^2-2Rr\cos(\theta)+R^2)^{\frac{3}{2}}}\sin(\theta) \space d\theta \space dr.$$
I hit a brick wall upon trying to evaluate the integral - ordinarily I would use a substitution in the single integral case but am unsure of how to do so for a double integral when the variables are all mixed up. Am I on the right track?
I'm not sure that this integral converges, given that the other two diverge, does this formula apply to point charges or only to continuous charge distributions?