Hopefully Lionel will expand his answer because it looks elegant but I can't wholly make it work. In the meanwhile, here's a ham fisted approach:
From the standard identity $\nabla \cdot(\vec A \times \vec B)= \vec{B}\cdot (\nabla \times \vec{A}) - \vec{A}\cdot (\nabla \times \vec{B})$ and by Faraday's law to expand $\nabla \times \vec{B}$ on the rhs we get:
$$\mu_0^{-1} \nabla \cdot \left(\vec A \times \vec B \right)= \mu^{-1} {\left|\vec{B} \right|}^2 - \vec{A}\cdot\vec{J} - \epsilon_0 \vec{A}\cdot\partial_t \vec{A}
\tag{1}$$
Expanding $\epsilon_0\,\nabla\cdot(\phi\,\vec{E}) = -\epsilon_0 \nabla\phi\cdot{E} + \epsilon_0 \phi \nabla\cdot\vec{E}$ and using Gauss's law for electricity as well as $\vec{E} = -\nabla \phi - \partial_t \vec{E}$ we get:
$$\epsilon_0 \nabla\cdot \left(\phi \vec{E} \right) = -\epsilon_0 \left|\vec{E}\right|^2 + \phi\,\rho + \epsilon_0 \vec{E}\cdot\partial_t\vec{A} \tag{2}$$
By applying the divergence theorem to a big spherical volume $V$ enclosing all of the current and charge to the left-hand side of both $\left(1\right)$ and $\left(2\right)$ you'll get surface integrals $\oint_{\partial V}\hat{n} \cdot \left(\vec{A}\times\vec{B} \right)\, {\rm d}S$ and $\oint_{\partial V}\hat{n}\cdot \left( \phi \vec{E} \right)\, {\rm d}S$. As in the comments, you need to assume decay behavior of $\vec{E}$, $\phi$, $\vec{B}$ and $\vec{A}$ their asymptotic dependencies of their magnitudes on the radius $R$ of the bign sphere for static conditions are at most $R^{-1}$ for $\phi$ and $\vec{A}$ and at most $R^{-2}$ for $\vec{E}$ and $\vec{B}$ (in general with dynamic problems with radiation, they all dwindle like $R^{-1}$). So in the static case, the surface integrals will vary like $R^{-1}$, hence vanish as the sphere grows big enough to enclose all space. On dropping the time varying terms on the right-hand side of $\left(1\right)$ and $\left(2\right)$ (for static conditions), you'll get the result you need.
As in the comments, there are easier – and more general – ways to study energy flows in the electromagnetic field. In the middle of my answer to the Physics SE question "How can Magnets be used to pick up pieces of metal when the force from a magnetic field does no work?" I show the standard method as in e.g. Griffiths or the literary horror of the symbol jungles of Jackson, "Classical Electrodynamics" (be sure to take a keenly sharpened machete with you to cut through all the equations).
The very best exposition I know of is in Chapter 27 of the second volume of the Feynman Lectures on Physics: the chapter called "Field Energy and Field Momentum". Feynman is every bit as mathematically rigorous as Jackson (probably better) and he studies the physics most carefully: discussing the often glossed over topics of the ambiguity in defining the field energy and flow as well as the locality of energy flows. The method in your question is interesting, though: I've not seen it before and it lets one split apart the $\frac{\epsilon}{2} \left| \vec{E}^2 \right|$ and $\frac{\mu}{2} \left|\vec{H}^2 \right|$ bits of the energy density to study the Lorentz-covariant $\left| \vec{E} \right|^2 - c^2 \left|\vec{B} \right|^2$.
Hope this all helps.
