Thanks for the update. "Irrotational" means curl = zero, so we compute the curl, set it equal to zero, and solve for the constants.
[tex]\vec f(x,y,z) = (axy+z^3) \, \vec\imath + (3x^2-z) \, \vec\jmath + (bz^2-y) \, \vec k[/tex]
[tex]\nabla \times \vec f(x,y,z) = \left(\dfrac{\partial(bz^2-y)}{\partial y} - \dfrac{\partial(3x^2-z)}{\partial z}\right) \, \vec\imath - \left(\dfrac{\partial(bz^2-y)}{\partial x} - \dfrac{\partial(axy+z^3)}{\partial z}\right) \, \vec\jmath \\ ~~~~~~~~~~~~ + \left(\dfrac{\partial(3x^2-z)}{\partial x} - \dfrac{\partial(axy+z^3)}{\partial y}\right) \, \vec k[/tex]
[tex]\nabla \times \vec f(x,y,z) = 3z^2 \, \vec\jmath + (6-a)x \, \vec k[/tex]
At this point, all we can conclude is that a = 6 to make the k-component vanish. Unfortunately there's no zeroing out the j-component, so the field as given is *not* irrotational.
As I mentioned in comments, if the i-component had instead been [tex]axy+bz^3[/tex], we would have ended up with
[tex]\nabla \times \vec f(x,y,z) = 3bz^2 \, \vec\jmath + (6-a)x \, \vec k[/tex]
in which case we would have b = 0.
I'll continue working with this "fixed" field to find Φ, if only to give you an idea of how to proceed. We want a scalar function Φ(x, y, z) such that the given vector field is the gradient of f. This would entail solving the partial differential equations,
[tex]\dfrac{\partial\Phi}{\partial x} = axy+bz^3 = 6xy[/tex]
[tex]\dfrac{\partial\Phi}{\partial y} = 3x^2 - z[/tex]
[tex]\dfrac{\partial\Phi}{\partial z} = bz^2-y = -y[/tex]
Let's start with the last equation. Integrating both sides with respect to z yields
[tex]\Phi(x,y,z) = \displaystyle \int (-y) \, dz = -yz + g(x,y)[/tex]
Now differentiate both sides with respect to y :
[tex]\dfrac{\partial\Phi}{\partial y} = -z + \dfrac{\partial g}{\partial y} = 3x^2 - z \implies \dfrac{\partial g}{\partial y} = 3x^2[/tex]
Integrate both sides of the latter PDE with respect to y to solve for g :
[tex]g(x,y) = \displaystyle \int 3x^2 \, dy = 3x^2y + h(x)[/tex]
Now we differentiate Φ with respect to x :
[tex]\Phi(x,y,z) = -yz + 3x^2y + h(x) \\\\ \dfrac{\partial\Phi}{\partial x} = 6xy + \dfrac{dh}{dx} = 6xy \implies \dfrac{dh}{dx} = 0[/tex]
Integrate to solve for h :
[tex]h(x) = \displaystyle \int 0 \, dx = C[/tex]
where C is an arbitrary constant.
So the scalar function whose gradient is our "fixed" field f is
[tex]\Phi(x,y,z) = -yz + 3x^2y + C[/tex]
