En las páginas 498-499 del libro de Lee "Introduction to Smooth Manifolds" se puede encontrar una discusión de una técnica para encontrar variedades integrales basada en la demostración del teorema de Frobenius. La idea básica es considerar una proyección $\pi \colon \mathbb{R}^3 \rightarrow \mathbb{R}^2$ tal que $d\pi|_{\mathscr{D}}$ es un isomorfismo cerca de $p$ y encontrar campos vectoriales $V,W$ que abarcan $\mathscr{D}$ cerca de $p$ que son $\pi$ -relacionados con los campos vectoriales de coordenadas estándar $\partial_x,\partial_y$ en $\mathbb{R}^2$ . Desde $\partial_x,\partial_y$ desplazarse y $d\pi|_{\mathscr{D}}$ es un isomorfismo cerca de $p$ los campos vectoriales $V,W$ también conmutan y podemos utilizar sus flujos para construir variedades integrales de $\mathscr{D}$ .
En su caso, podemos trabajar con la proyección sobre el $xy$ y sustituir los campos vectoriales $X,Y$ con los campos vectoriales
$$ 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. $$
Claramente, $V$ y $W$ son linealmente independientes y también abarcan $\mathscr{D}$ pero como $\pi_{*}(V) = \partial_x, \pi_{*}(W) = \partial_y$ también tenemos que $V,W$ conmutan (a diferencia de $X,Y$ para lo cual $[X,Y] = Y$ ).
El flujo de $V$ viene dado por
$$ \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) $$
mientras que el flujo de $W$ viene dado por
$$ \varphi^W_t(x_0,y_0,z_0) = \left( x_0, y_0 + t, \frac{z_0 x_0}{x_0 + t} \right). $$
Por lo tanto, el mapa $\Phi \colon \mathbb{R}^3 \rightarrow \mathbb{R}^3$ dada por
$$ \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) $$
nos dará un gráfico local de $\mathbb{R}^3$ con $\Phi(0,0,0) = (1,1,1)$ y tal que $(u,v) \mapsto \Phi(u,v,0)$ es una colector integral de $\mathscr{D}$ de paso $(1,1,1)$ . Para describir implícitamente este colector, basta con resolver para $(u,v,w)$ en términos de $(x,y,z)$ y luego $w = 0$ nos dará el colector integral. En este caso, obtenemos
$$ w = zx(x + y - 1) - 1 = 0. $$