General solution for the system of PDEs from the curl of a vector field equaling another In my vector calculus class, when we were introduced to the curl operator the professor gave us this example:
Is it possible to find a vector field $\mathbf{G}$ such that 
$$\mathbf{F} = \nabla \times {\mathbf{G}}?$$
As a motivation for the identity that the divergence of a curl of a vector field is $0$. 
My question is:

Given that the identity holds true, how do you solve for a general solution for the component functions of $\mathbf{G}$, given $\mathbf{F}$?

The best I can do is find one or two solutions. For example, if $\mathbf{F}$ = $\langle-y,-z,-x\rangle$, a solution I find is $\mathbf{G}$ = $\langle xy,0,-\frac{1}{2}y^2+xz\rangle$. I had to go through a system of partial differential equations and made it work, but can the general solution be written explicitly?
Thanks!

$\newcommand{\F}{\mathbf{F}}\newcommand{\p}{\mathbf{p}}$If we consider the unbounded $\mathbb{R}^3$ case, there is a path integral formula to construct the right inverse of the curl operator for divergence free vector field 
$$
\mathcal{R}(\F) = -(\p - \p_0)\times \int^1_0 \F\Big(\p_0 + t(\p- \p_0)\Big)t\,dt.\tag{1}
$$
A proof of this formula can be found here, for more discussion the author pointed to Spivak's book Calculus on Manifolds.
In your case, letting $\p_0 = \langle 0,0,0\rangle $, it is 
$$
\mathcal{R}(\F) = \langle x,y,z\rangle \times \int^1_0 \langle y,z,x\rangle t^2 dt
= -\frac{1}{3}\langle z^2-xy,x^2-yz,y^2-xz\rangle .$$
You can check that this is indeed the right inverse
$$
\nabla \times \mathcal{R}(\F) = \F = -\langle y, z, x\rangle.
$$
Notice the vector field above is difference from yours, because essentially $\mathcal{R}(\F) + \nabla \phi$ is also the answer for smooth $\phi$.

Further checking: the difference between the right inverse $\mathcal{R}(\F) $ above and your potential field $\langle xy,0,−y^2/2+xz\rangle$ is
$$
\mathbf{A} = \langle \frac{2}{3}x y +\frac{1}{3}z^2, \frac{1}{3}(x^2-y z), -\frac{1}{6}y^2+\frac{2}{3}xz\rangle = \nabla \left(\frac{1}{3} x^2y + \frac{1}{3}xz^2 - \frac{1}{6}y^2 z\right),
$$
indeed a gradient. 
If you wanna eliminate the huge kernel of curl operator, what you can do is (1) choosing a gauge (pinning down the divergence of $\mathbf{G}$), and/or (2) specifying boundary condition restricted on a simply-connected domain $\Omega$ of $\mathbb{R}^3$, by posing the following boundary value problem on some $\Omega$:
$$\left\{
\begin{aligned}
\nabla \times \mathbf{G} &= \F\quad \text{ in }\Omega,
\\
\nabla \cdot \mathbf{G} &= g \quad \text{ in }\Omega,
\\
\mathbf{G} \cdot \mathbf{n} &= 0 \quad \text{ on }\Gamma.
\end{aligned}
\right.$$
You can check the construction of yours $\langle xy,0,−y^2/2+xz \rangle $ is not divergence free, while formula (1) automatically gives you the divergence free potential field. For unbounded $\mathbb{R}^3$, choosing a gauge $g$ will get you the result you want.
