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.
