Solution (source code)

= Solution

Integrating $\psi_x=-(2xy+z)$ with respect to $x$ gives $\psi=-x^2y-xz+h(y,z)$. Matching $\psi_y=-(x^2+2y)$ requires $h=-y^2+j(z)$, and matching $\psi_z=-x$ requires $j$ constant. A <scalar potential> with the requested sign is
$$
\boxed{\psi=-x^2y-xz-y^2+C.}
$$
Thus the <force field> is a <conservative vector field> on all of $\mathbb R^3$.