The example is interesting because $H$ is not bounded on $\mathbb R^3$
and the usual integral formulas from the Helmholtz decomposition cannot be applied.
A solenoidal field $G$ is not one that has grad$\,G=0$ as you wrote but div$\,G=0\,.$ From
$$
\operatorname{div}H=2(xy+yz+zx)\,,\quad\operatorname{curl}H=-\left(\begin{smallmatrix}y^2\\z^2\\x^2\end{smallmatrix}\right)
$$
it follows that you are looking for fields $G$ and $F$ that have
\begin{align}
\operatorname{div}F&=2(xy+yz+zx)\,,\quad \operatorname{curl}G=-\left(\begin{smallmatrix}y^2\\z^2\\x^2\end{smallmatrix}\right)\,.
\end{align}
Since $\operatorname{curl}F=0$ it is of the form $F=\operatorname{grad}\Phi$
with some scalar function $\Phi\,.$ What you have to solve now is the Poisson equation
\begin{align}
\operatorname{div\,grad}\Phi=\Delta\Phi&=2(xy+yz+zx)\,.
\end{align}
Once this is done it follows that $G:=H-\operatorname{grad}\Phi$ is divergence free.
A simpler approach is the following:
Using the anti-curl formula from this post yields a field $G'$
\begin{align}
G'&:=\left[-\int_0^1t\left(\begin{smallmatrix}(ty)^2\\(tz)^2\\(tx)^2 \end{smallmatrix}\right)\,dt\right]\times\left(\begin{smallmatrix}x\\y\\z\end{smallmatrix}\right)=-\frac{1}{4}\begin{pmatrix}y^2\\z^2\\x^2\end{pmatrix}\times\begin{pmatrix}x\\y\\z\end{pmatrix}=-\frac{1}{4}\begin{pmatrix}z^3-x^2y\\x^3-y^2z\\y^3-z^2x\end{pmatrix}\\
&=-\frac{1}{4}\begin{pmatrix}z^3\\x^3\\y^3\end{pmatrix}+\frac{1}{4}H
\end{align}
which has the desired $\operatorname{curl}G'=-\left(\begin{smallmatrix}y^2\\z^2\\x^2\end{smallmatrix}\right)$ but is not divergence free: ($\operatorname{div}G'=(xy+yz+zx)/2$).
Still, the ansatz
$$
F'=H-G'=\frac{3}{4}H+\frac{1}{4}\begin{pmatrix}z^3\\x^3\\y^3\end{pmatrix}=\frac{1}{4}\begin{pmatrix}z^3+3x^2y\\x^3+3y^2z\\y^3+3z^2x\end{pmatrix}
$$
leads to a field $F'$ that has
$$
\operatorname{div}F'=\frac{3}{2}(xy+yz+zx)\,.
$$
Therefore,
$$
\boxed{\quad
F:=\frac{1}{3}\begin{pmatrix}z^3+3x^2y\\x^3+3y^2z\\y^3+3z^2x\end{pmatrix}
=\frac{1}{3}\begin{pmatrix}z^3\\x^3\\y^3\end{pmatrix}+H\quad}
$$
has the desired divergence $\operatorname{div}F=2(xy+yz+zx)\,.$
It can be checked that $F$ is rotation free.
Finally, $G$ is
$$\boxed{\quad
G=H-F=\frac{1}{3}\begin{pmatrix}z^3\\x^3\\y^3\end{pmatrix}\quad}
$$
which is obviously divergence-free.
The function $\Phi$ whose gradient is $F$ is
$$\Phi(x,y,z)=\frac{x^3y+y^3z+z^3x}{3}\,.$$