Chapter II, problem 2.8(b)  of Zachmanoglu and Thoe
Find the general solution to 
dx/( xy)= dy/(-xz)=(dz/((xy((x^2)+(y^2)))))
In class, we found two first integrals, and the region Ω where the vector field V is nonzero.  Indeed, Ω is (R^3) with the three coordinate axes removed.
In class I found two first integrals
u1(x,y,z) = x^2 + y^2 and
u2(x,y,z) = (x^2 + y^2) x^2  - z^2.

We can define a third First integral, u3 = u1^2-u2:
u3(x,y,z) = (x^2 + y^2) y^2 + z^2.

u1[x_, y_, z_] := x^2 + y^2;
u2[x_, y_, z_] :=(x^2 + y^2)x^2 - z^2;
u3[x_,y_,z_] := (x^2+y^2)y^2 + z^2;  (* u3 = u1^2-u2 *)
{x0, y0, z0} = {0,-1,1};
c1 = u1[x0, y0, z0];
c2 = u2[x0, y0, z0];
c3 = u3[x0, y0, z0];
show1 = True;
show2 = True;
show3 = True;
p1 = ContourPlot3D[u1[x,y,z] == c1, {x, -2,2}, {y, -2,2}, {z, -2,2}, ContourStyle -> Directive[Orange,Opacity[.7],Specularity[White,20]], AxesLabel ->{x,y,z}, Mesh -> None];p2 = ContourPlot3D[u2[x,y,z] == c2, {x, -2,2}, {y, -2,2}, {z, -2,2}, ContourStyle -> Directive[Darker[Green],Opacity[.7],Specularity[White,20]], AxesLabel ->{x,y,z}, Mesh -> None];p3 = ContourPlot3D[u3[x,y,z] == c3, {x, -2,2}, {y, -2,2}, {z, -2,2}, ContourStyle -> Directive[Cyan,Opacity[.7],Specularity[White,20]], AxesLabel ->{x,y,z}, Mesh -> None];ps = Graphics3D[Sphere[{x0, y0, z0}, .1]] ;
figs = {ps};
If[show1, AppendTo[figs, p1]];
If[show2, AppendTo[figs, p2]];
If[show3, AppendTo[figs, p3]];
Show[ figs, PlotRange ->{{-2,2},{-2,2}, {-2,2}}, Axes -> True,AxesLabel ->{x,y,z}]

v[x_, y_, z_] := {y z, -x z, x y(x^2+y^2)};gradu1[x_,y_, z_] := {2x, 2y, 0};
gradu2[x_,y_,z_] :={4 x^3 + 2x y^2, 2 x^2 y, -2z}
gradu3[x_,y_,z_] :={2 x y^2, 2 x^2 y + 4 y^3, 2z}

Simplify[v[x,y,z]. gradu1[x,y,z] ]
Simplify[v[x,y,z]. gradu2[x,y,z] ]
Simplify[v[x,y,z]. gradu3[x,y,z] ]

Simplify@Cross[gradu1[x,y, z], gradu2[x,y, z]]
Simplify@Cross[gradu1[x,y, z], gradu3[x,y, z]]
Simplify@Cross[gradu2[x,y, z], gradu3[x,y, z]]

v[x,y,z]

For each pair the cross product is nonzero precisely where v is nonzero.   Thus each pair is a set of functionally independent first integrals on Omega = R^3 with the three coordinate axes removed.
Note that grad u2 x grad u3 is a scalar valued function, 8(x^2 + y^2), times v. In the other two cases that scalar function is a constant (4 or -4), which is not to be expected in general.

Note that there is a nontrivial function mapping R^3 to R, namely F(a,b,c) = a^2 - b - c, such that F(u1(x,y,z), u2(x,y,z), u3(x,y,z))  = 0 for all (x,y,z) in Omega.