Electrostatic Potential Energy Integral 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? 
 A: From Griffith section 2.4.4 comments on Electrostatic Energy, you can get your answer. 
If you consider point charges, then actually, this integral is related with self-energy which is infinite at usual, 
to make finite we often introduce cutoff radius $\delta$.  (In particle physics, we often use bare and renormalized terminology, renormalization is a some process make infinte to finite) 
The relevant integral is well describe in Griner's Electrodynamics and Jackson's ch1. 
These two textbook contains both calculation and its physical interpretation as well. 
A: The Poynting formula for electrostatic energy in volume $V$
$$
E = \int_V \frac{1}{2}\epsilon_0 E^2 dV
$$
can be derived from the Coulomb law only for cases where the field acting on the particles is defined everywhere. However, point particle has infinite charge density at the point it is present and the field is not defined at that point. So the derivation fails.
For two point particles at rest, the work necessary to bring these particles to their positions $\mathbf r_1,\mathbf r_2$ is known to be
$$
W = \frac{1}{4\pi\epsilon_0}\frac{q_1q_2}{|\mathbf r_1- \mathbf r_2|}
$$
If you want to express this energy in terms of EM fields only, this can be written as
$$
\int_{whole~space} \epsilon_0\mathbf E_1(\mathbf x) \cdot \mathbf E_2(\mathbf x) \,d^3\mathbf x
$$
where $\mathbf E_1(\mathbf x) = -\nabla \phi_1(\mathbf x)$ is field due to the first particle
and $\mathbf E_2(\mathbf x)=-\nabla \phi_2$ is field due to the second particle.
Proof:
The potential $\phi_1$ is
$$
\mathbf \phi_1(\mathbf x) = \frac{1}{4\pi\epsilon_0}\frac{q_1}{|\mathbf x - \mathbf r_1|}
$$
and the potential $\phi_2(\mathbf x)$ is
$$
\mathbf \phi_2(\mathbf x) = \frac{1}{4\pi\epsilon_0}\frac{q_2}{|\mathbf x - \mathbf r_2|}.
$$
The integral becomes
$$
\int_{whole~space} \epsilon_0\mathbf E_1(\mathbf x) \cdot \mathbf E_2(\mathbf x) \,d^3\mathbf x = \int_{whole~space} \epsilon_0\nabla\phi_1(\mathbf x) \cdot \nabla \phi_2(\mathbf x) \,d^3\mathbf x = 
$$
$$
= \int_{whole~space} \epsilon_0\nabla\cdot( \phi_1 \nabla \mathbf \phi_2 )\,d^3\mathbf x -\int_{whole~space} \epsilon_0\phi_1 \Delta \phi_2\,d^3\mathbf x.
$$
For electrostatic field, the first integral is zero (this can be shown using the Gauss theorem). For the second potential, the Poisson equation 
$$
\Delta \phi_2 = -\frac{q_2}{\epsilon_0}\delta(\mathbf x - \mathbf r_2)
$$
holds so we arrive at the integral
$$
\int_{whole~space} \frac{1}{4\pi\epsilon_0}\frac{q_1}{|\mathbf x - \mathbf r_1|}\frac{q_2}{\epsilon_0}\delta(\mathbf x - \mathbf r_2)\,d^3\mathbf x
$$
which has the value
$$
\frac{1}{4\pi\epsilon_0}\frac{q_1q_2}{|\mathbf r_2 - \mathbf r_1|},
$$
which is the same as $W$ above.
This formula for EM energy has general version for time-dependent fields
$$
E_{em} = \int \epsilon_0\mathbf E_1\cdot\mathbf E_2 + \frac{1}{\mu_0}\mathbf B_1\cdot \mathbf B_2\,d^3\mathbf x
$$
In case more particles are involved, similar formulae can be derived, with summation over each pair of particles.
R. C. Stabler, A Possible Modification of Classical Electrodynamics, Physics Let-
ters, 8, 3, (1964), p. 185-187.
http://dx.doi.org/10.1016/S0031-9163(64)91989-4
J. A. Wheeler, R. P. Feynman, Classical Electrodynamics in Terms of Direct
Interparticle Interaction, Rev. Mod. Phys., 21, 3, (1949), p. 425-433.
http://dx.doi.org/10.1103/RevModPhys.21.425
J. Frenkel, Zur Elektrodynamik punktförmiger Elektronen, Zeits. f. Phys., 32, (1925), p. 518-534.
http://dx.doi.org/10.1007/BF01331692
