Your approach is absolutely correct, but you're not using the formula of energy density correctly. Energy per unit volume stored in electric field at some point $P$ is: \begin{equation}\frac{\epsilon}{2}E_p^2 \end{equation} where $E_p$ is the electric at $P$ and $\epsilon$ is the permittivity of the material present at $P$.

That being said, when you apply this formula to find the energy stored inside the dielectric, the field $E_p$ will be the field inside the dielectric($E_d$) and not the field in the vacuum($E$). \begin{equation} E_d=\frac{\epsilon_0E}{\epsilon}\quad ...(\star)\quad(\epsilon:\text{permittivity of dielectric}) \end{equation}Now, I'll do the same mathematics as you did but I'll use the formula correctly.

I'll make the assumptions you made: The dielectric is small enough the electric field can be assumed to be constant over the entire volume($V$) of the dielectric.

Electrostatic potential energy before moving the dielectric $W$ (includes the energy of dielectric at $P$and energy of the space where the dielectric will be moved): \begin{equation} \frac{W}{V}=\frac{\epsilon}{2}E_d^2+\frac{\epsilon}{2}(E^2+2\vec{E}.d\vec{E})= \frac{\epsilon_0^2}{2\epsilon}E^2 + \frac{\epsilon_0}{2}E^2 +\frac{\epsilon_0}{2}dE^2 \end{equation} $\quad$ (using ($\star$). Also $\vec{E}.d\vec{E}=\frac{1}{2}d(\vec{E}.\vec{E})=\frac{1}{2}dE^2$ ) $\quad$ where $d\vec{E}$ is the change in electric field when we move $d\vec{r}$ from $P$. Now, potential energy after moving the dielectric by $d\vec{r}$(includes the potential energy stored inside dielectric after it has been moved and the energy stored in space where the dielectric was initially): \begin{equation} \frac{W+dW}{V}=\frac{\epsilon_0}{2}E^2+\frac{\epsilon}{2}(E_d^2+2\vec{E_d}.d\vec{E_d})=\frac{\epsilon_0}{2}E^2+\frac{\epsilon_0^2}{2\epsilon}E^2+\frac{\epsilon_0 ^2}{2 \epsilon}dE^2 \end{equation} Finally, \begin{equation} \frac{dW}{V}=\frac{\epsilon_0 ^2}{2 \epsilon}dE^2-\frac{\epsilon_0}{2}dE^2=\frac{\epsilon_0}{2}\left(\frac{\epsilon_0}{\epsilon}-1 \right)dE^2 \Rightarrow dW=\frac{\epsilon_0 V}{2}\left(\frac{\epsilon_0}{\epsilon}-1 \right)dE^2 \end{equation} This change in potential energy is equal to the negative of the work done by the conservative force (Electric force, $\vec{F}_e$ here): \begin{equation} -\vec{F}_e.d\vec{r}=\frac{\epsilon_0 V}{2}\left(\frac{\epsilon_0}{\epsilon}-1 \right)dE^2 \Rightarrow \vec{F}_e.d\vec{r}=\frac{\epsilon_0 V}{2}\left(1- \frac{\epsilon_0}{\epsilon} \right)dE^2= (\Lambda).dE^2\quad ...(\star\star) \end{equation} $\quad$ where $\Lambda$ is some positive real number.

Now if $d\vec{r}$ corresponds to a direction along which $E$ increases, RHS of the $(\star\star)$ is positive because $d(E^2)$ is positive. In which case, LHS also needs to be positive for the equation to hold $\therefore$ Electrostatic force will be along the direction of $d\vec{r}$, the direction of increase in $E$, for the dot product to stay positive.

You can also write: $\vec{F}_e=\Lambda \vec{\nabla}E^2$ and deduce the same result, though more precisely.

Your approach is absolutely correct, but you're not using the formula of energy density correctly. Energy per unit volume stored in electric field at some point $P$ is: \begin{equation}\frac{\epsilon}{2}E_p^2 \end{equation} where $E_p$ is the electric field at $P$ and $\epsilon$ is the permittivity of the material present at $P$.

That being said, when you apply this formula to find the energy stored inside the dielectric, the field $E_p$ will be the field inside the dielectric($E_d$) and not the field in the vacuum($E$). \begin{equation} E_d=\frac{\epsilon_0E}{\epsilon}\quad ...(\star)\quad(\epsilon:\text{permittivity of dielectric}) \end{equation}Now, I'll do the same mathematics as you did but I'll use the formula correctly.

I'll make the assumptions you made: The dielectric is small enough, therefore electric field can be assumed to be constant over the entire volume($V$) of the dielectric.

Electrostatic potential energy before moving the dielectric ($W$ is the energy, I divided it by Volume of the dielectric. We can rearrange the terms later) (includes the energy of dielectric at $P$and energy of the space where the dielectric will be moved): \begin{equation} \frac{W}{V}=\frac{\epsilon}{2}E_d^2+\frac{\epsilon}{2}(E^2+2\vec{E}.d\vec{E})= \frac{\epsilon_0^2}{2\epsilon}E^2 + \frac{\epsilon_0}{2}E^2 +\frac{\epsilon_0}{2}dE^2 \end{equation} $\quad$ (using ($\star$). $\quad$ Also $\vec{E}.d\vec{E}=\frac{1}{2}d(\vec{E}.\vec{E})=\frac{1}{2}dE^2$ ) $\quad$ where $d\vec{E}$ is the change in electric field when we move $d\vec{r}$ away from $P$.

Now, potential energy after moving the dielectric by $d\vec{r}$(includes the potential energy stored inside dielectric after it has been moved and the energy stored in space where the dielectric was initially): \begin{equation} \frac{W+dW}{V}=\frac{\epsilon_0}{2}E^2+\frac{\epsilon}{2}(E_d^2+2\vec{E_d}.d\vec{E_d})=\frac{\epsilon_0}{2}E^2+\frac{\epsilon_0^2}{2\epsilon}E^2+\frac{\epsilon_0 ^2}{2 \epsilon}dE^2 \end{equation} Finally, \begin{equation} \frac{dW}{V}=\frac{\epsilon_0 ^2}{2 \epsilon}dE^2-\frac{\epsilon_0}{2}dE^2=\frac{\epsilon_0}{2}\left(\frac{\epsilon_0}{\epsilon}-1 \right)dE^2 \Rightarrow dW=\frac{\epsilon_0 V}{2}\left(\frac{\epsilon_0}{\epsilon}-1 \right)dE^2 \end{equation} This change in potential energy is equal to the negative of the work done by the conservative force (Electric force, $\vec{F}_e$ here): \begin{equation} -\vec{F}_e.d\vec{r}=\frac{\epsilon_0 V}{2}\left(\frac{\epsilon_0}{\epsilon}-1 \right)dE^2 \Rightarrow \vec{F}_e.d\vec{r}=\frac{\epsilon_0 V}{2}\left(1- \frac{\epsilon_0}{\epsilon} \right)dE^2= (\Lambda).dE^2\quad ...(\star\star) \end{equation} $\quad$ where $\Lambda$ is some positive real number.

Now if $d\vec{r}$ corresponds to a direction along which $E$ increases, RHS of the $(\star\star)$ is positive because $d(E^2)$ is positive. In which case, LHS also needs to be positive for the equation to hold $\therefore$ Electrostatic force will be along the direction of $d\vec{r}$, the direction of increase in $E$, for the dot product to stay positive.

You can also write: $\vec{F}_e=\Lambda \vec{\nabla}E^2$ and deduce the same result, though more precisely.