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. $$