Solution (source code)

= Solution

Integrating the normal Maxwell traction (the radial parts cancel) gives
$$
F_z=\frac{q^2}{16\pi\epsilon_0d^2}\frac{\epsilon_r-1}{\epsilon_r+1}.
$$
It points toward the dielectric when $\epsilon_r>1$. This equals the Coulomb force from the image charge, as it must.