Let $X,Y \in \mathfrak{X}(\mathbb{R}^3)$ be vector fields on $\mathbb{R}^3$ given by $$X = x \displaystyle{\frac{\partial}{\partial x}} + \displaystyle{\frac{\partial}{\partial y}} + x(y + 1) \displaystyle{\frac{\partial}{\partial z}}, \ Y = \displaystyle{\frac{\partial}{\partial x}} + y\displaystyle{\frac{\partial}{\partial z}}.$$
How can I determine $U,V \in \mathfrak{X}(\mathbb{R}^3)$ such that $[U,V] = 0$ and the distribution $D(U,V)$ generated by $U,V$ coincides with the distribution $D(X,Y)$ generated by $X,Y$?
I know that $D(X,Y)$ is integrable (as $[X,Y] = -Y$ we have $D(X,Y)$ involutive, hence Frobenius Theorem ensures $D(X,Y)$ integrable), but I do not know if it helps. I've tried to deduce conditions on $U$ and $V$ to have both $[U,V] = 0$ and $D(X,Y) = D(U,V)$ and it led me nowhere.
Any ideas on how to tackle such a problem?
As always, I find it easiest to work with the $1$-form(s) annihilating the distribution. In your case, you can check easily that $$\omega = y\,dx + x\,dy - dz$$ satisfies $\omega(X)=\omega(Y)=0$. Certainly $\omega=0$ is integrable since $\omega = d(xy-z)$. (This tells you that the integral manifolds are the surfaces $xy-z= c$ for different constants $c$.)
Can you now easily find the $U,V$ you want?