Finding an explicit integral manifold Let $\mathscr{D}$ be the $2$-dimensional distribution on $M:=\{(x,y,z)\in\mathbb{R}^3:x,y,z>0\}$ generated by $\{X,Y\}$ where $X=x\partial_x - 2y\partial_y, Y=xy\partial_y - xz\partial _z$.
By a direct computation, I showed that $\mathscr{D}$ is involutive. Hence, by Frobenius theorem, for every point $p$ of $M$, there exists an integral manifold of $\mathscr{D}$ passing through $p$.
Let's pick $(1,1,1)\in M$.
How do I construct an explicit integral manifold of $\mathscr{D}$ passing through $(1,1,1)$?
 A: You can find a discussion of a technique for finding integral manifolds based on the proof of the Frobenius theorem in pages 498-499 of Lee's "Introduction to Smooth Manifolds" book. The basic idea is to consider a projection $\pi \colon \mathbb{R}^3 \rightarrow \mathbb{R}^2$ such that $d\pi|_{\mathscr{D}}$ is an isomorphism near $p$ and find vector fields $V,W$ that span $\mathscr{D}$ near $p$ which are $\pi$-related to the standard coordinate vector fields $\partial_x,\partial_y$ on $\mathbb{R}^2$. Since $\partial_x,\partial_y$ commute and $d\pi|_{\mathscr{D}}$ is an isomorphism near $p$, the vector fields $V,W$ also commute and we can use their flows to construct integral manifolds of $\mathscr{D}$.
In your case, we can work with the projection on the $xy$ plane and replace the vector fields $X,Y$ with the vector fields
$$ V = \frac{1}{x} X + \frac{2}{x^2} Y = \partial_x - \frac{2z}{x} \partial_y, \\
W = \frac{1}{xy} Y = \partial_y - \frac{z}{y} \partial_z. $$
Clearly, $V$ and $W$ are linearly independent and also span $\mathscr{D}$ but since $\pi_{*}(V) = \partial_x, \pi_{*}(W) = \partial_y$ we also have that $V,W$ commute (unlike $X,Y$ for which $[X,Y] = Y$).
The flow of $V$ is given by
$$ \varphi^V_t(x_0,y_0,z_0) = \left( x_0 + t, y_0, \frac{z_0 x_0^2}{(x_0 + t)^2} \right) $$
while the the flow of $W$ is given by
$$ \varphi^W_t(x_0,y_0,z_0) = \left( x_0, y_0 + t, \frac{z_0 x_0}{x_0 + t} \right). $$
Hence, the map $\Phi \colon \mathbb{R}^3 \rightarrow \mathbb{R}^3$ given by
$$ \Phi(u,v,w) = (\varphi^V_{u} \circ \varphi^W_{v})(1,1,1 + w) = \left( 1 + v, 1 + u, \frac{1 + w}{(1 + v)(1 + v + u)} \right) $$
will give us a local chart of $\mathbb{R}^3$ with $\Phi(0,0,0) = (1,1,1)$ and such that $(u,v) \mapsto \Phi(u,v,0)$ is an integral manifold of $\mathscr{D}$ passing through $(1,1,1)$. To describe this manifold implicitly, we just need to solve for $(u,v,w)$ in terms of $(x,y,z)$ and then $w = 0$ will give us the integral manifold. In this case, we get
$$ w = zx(x + y - 1) - 1 = 0. $$
A: To find an integral manifold of involutive distribution $\mathscr{D}$, we can construct the differential ideal, where differential forms annihilate $\mathscr{D}$. In your case, for dimensional reason, the ideal is just generated by a one form $$\omega=P\mathrm{dx}+Q\mathrm{dy}+R\mathrm{dz}$$
Because $\omega$ annihilate $\mathscr{D}$, $(P,Q,R)$ is parellel to the normal vector $$(x,-2y,0)\times(0,xy,-xz)=(2xyz,x^2z,x^2y)$$
and we can take $\omega=2xyz\mathrm{dx}+x^2z\mathrm{dy}+x^2y\mathrm{dz}$, which is an exact form $$\omega=d(x^2yz)$$
The integral manifold is given by level set $$x^2yz=c$$
where $c=1$ for passing through $(1,1,1)$.
