We can make use of the multivariable chain rule (see this answer) which states that given a function $h: \mathbb{R}^n \to \mathbb{R}^m$ and a function $f: \mathbb{R}^m \to \mathbb{R}$ we have
$$
\nabla (f \circ h) = \left(h'\right)^T \cdot \left(\nabla f \circ h\right)
$$
where "$\cdot$" represents matrix multiplication.
In our case, we see that if we define $h(x,y) = (xy+x^2,xy-y^2)$ then it follows from the definition of $g$ that
$$g(x,y) = (f \circ h)(x,y)= f(xy+x^2,xy-y^2) \tag{1}$$
which means that calculating $\nabla g(3,-3)$ is equal to calculating $\nabla (f \circ h)(3,-3) $.
Now, by definition, we know that
$$
h' = \begin{pmatrix}
\frac{\partial}{\partial x}(xy+x^2) & \frac{\partial}{\partial y}(xy+x^2)\\
\frac{\partial}{\partial x}(xy-y^2) & \frac{\partial}{\partial y}(xy-y^2)
\end{pmatrix} = \begin{pmatrix}
y+2x & x\\
y & x-2y
\end{pmatrix}
$$
which means that the transposed matrix $\left(h'\right)^T$ is $\begin{pmatrix}
y+2x & y\\
x & x-2y
\end{pmatrix}
$. If we then evaluate this matrix at $(3,-3)$ we get
$$
\left(h'\right)^T(3,-3) = \begin{pmatrix}
-3+2(3) & -3\\
3 & 3-2(-3)
\end{pmatrix} = \begin{pmatrix}
-3+6 & -3\\
3 & 3+6
\end{pmatrix} = \begin{pmatrix}
3 & -3\\
3 & 9
\end{pmatrix}
$$
On the other hand, we see that
$$
h(3,-3) = (3(-3)+3^2, 3(-3) -(-3)^2) = (-9+9, -9-9) = (0,-18)
$$
which tells us that
$$
\left(\nabla f \circ h\right)(3,-3) = \nabla f\left(h(3,-3)\right) = \nabla f\left(0,-18\right) = \begin{pmatrix} -2 \\3 \end{pmatrix}
$$
using the convention of the gradient as a column vector.
Finally, putting all this together tells us that
$$
\nabla g(3,-3) = \nabla (f \circ h)(3,-3) = \left[\left(h'\right)^T(3,-3)\right] \cdot \left[\left(\nabla f \circ h\right)(3,-3)\right] = \begin{pmatrix}
3 & -3\\
3 & 9
\end{pmatrix} \cdot \begin{pmatrix} -2 \\3 \end{pmatrix} = \begin{pmatrix} -2(3)+3(-3) \\-2(3)+3(9) \end{pmatrix}= \begin{pmatrix} -6-9 \\-6+27 \end{pmatrix}= \begin{pmatrix} -15 \\21 \end{pmatrix} = -15 \boldsymbol{\hat\imath} + 21 \boldsymbol{\hat\jmath}
$$