DSolve[] can solve the system all at once:
DSolve[{Grad[f[x, y, z], {x, y, z}] == {y z, x z, x y}},
f[x, y, z], {x, y, z}]
(* {{f[x, y, z] -> x y z + C[1]}} *)
If you wish to perform the iterative approach you give in your notebook, then DSolve[] can perform the single-integral steps:
DSolve[{D[f[x, y, z], x] == y z}, f[x, y, z], {x, y, z}]
(* {{f[x, y, z] -> x y z + C[1][y, z]}} *)
Here C[1] is your g; that is C[1][y, z] is g[y, z]. This is an advantage of DSolve[] over Integrate here: Since the variables {x, y, z} are specified, DSolve[] makes the "constant of integration" depend on y and z, just what is needed for the next step.
You can then differentiate and integrate, eliminating a variable at each iteration as in your "Manual Method." I figure you would appreciate working out the rest of the steps yourself, but if you run into difficulty, you can post a reply and I will try to help.
