Deriving energy in Dielectrics Note: I am working in the Lorentz-Heaviside system and all the integrals are over the whole space.
Definitions:
$$\vec E= \vec E_f+\vec E_b$$
$$\phi=\phi_f+\phi_b$$
$$\vec D=\vec E+\vec P$$
$$\rho=\rho_f+\rho_b$$
Here the subscripts f and b refer to free and bound charges, thus $\vec E_f$ is the field created by free charge when all the charges are at the desired (final) positions.
$\\$
The total energy of the system is $$W=\frac12\int\rho\phi\, dV=\frac12\int E^2\, dV$$
But we are only interested in the energy that is required to move free charge from infinity to some desired position, thus we need to subtract the  energy needed to create the system of the bound charge itself, namely $$\frac12\int \rho_b\phi_b\,dV$$
Since this is the energy of creation of bound charge system, where $\phi_b$ is the final potential caused by this system.
One may think to also subtract the term of $\rho_b\phi_f\equiv\rho_f\phi_b$ but this is not correct since this is the interaction energy of two systems and must be included. Thus we get
$$W=\frac12\int\rho\phi\, dV-\frac12\int \rho_b\phi_b\,dV$$
$$W=\frac12\int(\rho_f+\rho_b)(\phi_f+\phi_b)\, dV -\frac12\int \rho_b\phi_b\,dV$$
$$W=\frac12\int\rho_f\phi_f\, dV+\int\rho_f\phi_b\, dV\tag{1}$$
$$W=\frac12\int\rho_f(\phi+\phi_b)\, dV$$
$$W=\frac12\int\vec E\cdot\vec D\, dV+\frac12\int\rho_f\phi_b\, dV\tag{2}$$
$\\$
We can also derive this as follows:
Consider that initially when the free charge is at infinity, we bring $dq_f$, the work needed for this will be
$$dw_1=dq_f\int_{\phi_{b_1}}^{\phi_{b_2}}\phi_b$$
For the second free charge
$$dw_2=dq_f\int_{\phi_{b_2}}^{\phi_{b_3}}\phi_b+dq_f\int_{\phi_{f_1}}^{\phi_{f_2}}\phi_f$$
Where the second integral comes due to the first charge.
Following the pattern, we get the net work to be
$$W=\int\rho_f\phi_b\, dV+\frac12\int\rho_f\phi_f\, dV$$
Here both the potentials are the final potentials of the respective configurations. This is the same as equation 1.
Thus the question is why are we getting an extra integral in equation 2? What is going wrong in this derivation?
 A: There are three two things wrong with your derivation.
No. 1 (not too bad) is that you assume $\vec{D}$ and $\vec{P}$ to be given as gradients of potentials, which in general is not given. (As already commented by basics) But ok, limiting to a special case should still yield the desired result, so that's not it. Silly me.
No. 2 is that you are missing the assumption of a linear medium that is usually used in deriving the integral $\int \vec{D} \cdot \vec{E} \,\text{d}V$ for the energy. In a Lagrangian formalism (and you might do the same with energy) the dielectric displacement is the conjugate quantity to the electric field, i.e. given a Lagrangian density $\mathcal{L}$ we have
$$
D_i = \frac{\partial \mathcal{L}}{\partial E_i}\,.
$$
For a linear isotropic medium, with the free-field Lagrangian being
$$
\mathcal{L} = \frac{1}{2} \varepsilon_0 \vec{E}^2
$$
you will find that $\vec{D} = \varepsilon_0 \vec{E}$ and thus indeed
$$
\mathcal{L} = \frac{1}{2} \vec{D} \cdot \vec{E}\,.
$$
In a nonlinear model this is not generally true, for example for the electrostatic Heisenberg-Euler Lagrangian
$$
\mathcal{L} = \frac{1}{2} \varepsilon_0 \vec{E}^2 + \frac{1}{2} \varepsilon_0 \vec{E}^2 \frac{\vec{E}^2}{E_0^2}
$$
($E_0$ being some arbitrary constant field)
you will find that
$$
\vec{D} = \varepsilon_0 \vec{E} + 2 \varepsilon_0 \vec{E} \frac{\vec{E}^2}{E_0^2}
$$
and subsequently
$$
\mathcal{L} = \frac{1}{2} \vec{D} \cdot \vec{E} - \frac{1}{2} \varepsilon_0 \vec{E}^2 \frac{\vec{E}^2}{E_0^2}\,.
$$
No. 3, the main point, is that you are starting with a wrong formula for the field energy. The integral
$$
W = \frac{1}{2} \int \varepsilon_0 \vec{E}^2 \,\text{d}V
$$
is the energy of a collection of charges in vacuum, which you for example obtain by bringing each charge into position one after the other from infinity and calculating the work you do each time. You cannot use this for the energy in a continuous medium because when building up the charge distribution in the medium you also do some internal work to move the polarization charges that is not accounted for in the above formula for $W$. What you do instead is consider the infinitesimal change of energy when changing the density of free charges by $\delta \rho$,
$$
\delta W = \int \delta\rho \phi \,\text{d}V\,.
$$
Through $\delta \rho = \nabla \cdot \left(\delta \vec{D}\right)$ and partial integration you obtain
$$
\delta W = \int \vec{E} \cdot \delta \vec{D} \,\text{d}V\,
$$
and work your way onwards from there. For example, for a linear medium you can integrate it to obtain the formula you wish to derive.
I highly recommend you have a look at chapter 4.7 in Jackson's "Classical Electrodynamics" (the "Bible"). Panofsky&Phillips, Landau&Lifshitz or others will serve you equally well.
P.S.: I like S.I. units in electrodynamics. :)
A: 
The total energy of the system is
$$W=\frac12\int\rho\phi\, dV=\frac12\int E^2\, dV$$

This is not what we usually call the "energy" of the system. Or at least it is not what we are usually interested in when we talk about the "energy." Usually, for linear media, we write the energy $W_{correct}$ as:
$$
W_{correct}=\frac12\int\rho_f\phi\, dV=\frac{1}{2}\int \vec E\cdot \vec D\, dV
$$

$$W=\frac12\int\vec E\cdot\vec D\, dV+\frac12\int\rho_f\phi_b\, dV\tag{2}$$


Thus the question is why are we getting an extra integral in equation 2? What is going wrong in this derivation?

Clearly, the problem is that you are starting from the wrong definition of "total energy."
Using your definition for $W$ and my definition for $W_{correct}$, we see that
$$
W - W_{correct} = \frac{1}{2}\int\rho_f\phi_b\, dV\;,
$$
which is exactly the extra part in your Eq (2) that you are worried about.
A: I think you messed up a little the path to go from work done to build the system of charges in the dielectrics to the definition of energy of the system.
From work to energy
Qualitatively, we can think at system energy as the total amount of work you do to build your final configuration (to be fair, the opposite of this work). Let's inspect the way you build your configuration:

*

*you start taking free charges from the infinity. It's ok, this way you're setting the additive constant to your energy, that has no meaning in classical physics;


*then, you take the first free charge $q_1^f$ and placed it in space: it's the first charge you move, this operation is "free", assuming that the dielectric is initially not charged or polarized;


*once you have placed the first charge, you can assess its influence in the surrounding space through the electric field, given as the sum of the contribution of the free charge and the polarization of the dielectric
$\mathbf{e}_1(\mathbf{r}) = \dfrac{\mathbf{d}_1(\mathbf{r})- \mathbf{p}_1(\mathbf{r})}{\epsilon_0}$,
being $\mathbf{d}_1(\mathbf{r})$ the displacement field generated by free charge $q_1^f$ and $\mathbf{p}_1(\mathbf{r})$ the polarization field created in the dielectric "as a response" to the free charge.


*when you take the second free charge $q_2^f$, a force
$\mathbf{F}_{21}(\mathbf{r}) = q_2^f \mathbf{e}_1(\mathbf{r})$
is acting on it, whose expression is given by the product of:

*

*the free charge $q_2^f$

*the overall electrical field (taking into account both the effects of free and bound charges) produced by the presence of the free charge $q_1^f$ in the dielectric.

Imagine you're dragging free charges with your hand: this is the force you feel on the free charge; what happens to bound charges can be interpreted as a consequence of free charges in a dielectric, and it's taken into account in the whole electric field $\mathbf{e}$, that is the sum of the contributions of free and bound charges.
Expression of the electrostatic energy of a system in a dielectric
Having in mind what we said in the previous paragraph, you can repeat the derivation you did in you question, to get the right expression of the electrostatic energy of a system

*

*for discrete charges:
$E = \dfrac{1}{2} \displaystyle\sum_{i,j\ne i} q^f_i \phi_j(\mathbf{r}_i)$


*for continuous charges, in infinite space (no boundary contributions):
$E = \dfrac{1}{2} \displaystyle \int_V \rho^f(\mathbf{r}) \phi(\mathbf{r}) = (\text{integration by parts}) = \dfrac{1}{2} \displaystyle \int_V \mathbf{d}(\mathbf{r}) \cdot \mathbf{e}(\mathbf{r}) $
