Ja so ähnlich ;).
Vorhin konnte ich leider keine Rechnungen anfügen, daher habe ich nur grob den Weg skizziert.
zum Gradientenfeld:
Rotation berechnet sich gemäß
https://de.wikipedia.org/wiki/Rotation_eines_Vektorfeldes#Definition_der_Rotation_in_kartesischen_Koordinaten
Das gibt bei uns (für p=1)
$$ \begin{pmatrix} x-x\\y-y\\z-z \end{pmatrix}=\begin{pmatrix} 0\\0\\0 \end{pmatrix} $$,
passt also.
Zum Potential:
Die 3 entstehenden Gleichungen lauten:
$$ \frac { \partial\varphi }{ \partial x }=yz+2x\\\frac { \partial\varphi }{ \partial y }=xz-2y\\\frac { \partial\varphi }{ \partial z }=xy\\$$
man nimmt z.B die erste Gleichung zur Hand und integriert:
$$ \varphi(x,y,z) =yzx+x^2 +C(y,z) $$
Diese Teillösung setzt man nun in die zweite DGl links ein, um den y-Anteil der Konstante C(y,z)
zu bestimmen:
$$ \frac { \partial\varphi }{ \partial y }= xz-2y\\zx+C_{y}(y,z)=xz-2y\\C_{y}(y,z)=-2y\\C(y,z)=-y^2+D(z)\\\to \varphi (x,y,z)= yzx+x^2-y^2+D(z)\\$$
Nun geht man denselben Weg, um D(z) zu bestimmen, mithilfe der dritten DGL.
Glücklicherweise erfüllt unsere Teillösung diese DGL bereits, daher ist D(z)=D eine echte Konstante.
Lösung:
$$ \varphi (x,y,z)= yzx+x^2-y^2+D\\$$