Submanifold Tangent to Given Vector Fields This is part of another old qualifying exam problem, so any help is appreciated! We are given two vector fields $X$ and $Y$ on $\mathbb{R}^3$ described by 
$$X = \frac{\partial}{\partial y} + z\frac{\partial}{\partial x}\qquad Y = \frac{\partial}{\partial z} + y\frac{\partial}{\partial x}$$
and we are asked to give a parameterization $(x(s,t), y(s,t), z(s,t))$ for the (unique) surface passing therough the point $p = (1,0,0)$ and tangent to the vector fields $X$ and $Y$ at each point, and to then give an equation of the form $F(x,y,z) = 0$ for this surface. 
I was able to determine that $[X, Y] = 0$, and that $X$ and $Y$ were linearly independent at $p$, so I computed the flows $\theta_t$ of $X$ and $\phi_s$ of $Y$ to be
$$\begin{align*}
\theta_t(x,y,z) &= (xe^{zt}, y + t, z)\\
\phi_s(x,y,z) &= (xe^{ys}, y, z + s)
\end{align*}$$
which gives us 
$$\theta_t\circ \phi_s = \phi_s\circ \theta_t = (xe^{zt + ys + st}, y + t, z + s),$$ but it is at this point that I get stuck. To solve this, I have tried to use the following theorem from Lee's textbook:

Theroem 9.46: Let $M$ be a smooth $n$ manifold, and let $(V_1, \ldots, V_k)$ be a linearly independent $k$-tuple of smooth commuting vector fields on an open subset $W\subseteq M$. For each $p\in W$, there exists a smooth coordinate chart $(U, (s^i))$ centered at $p$ such that $V_i = \frac{\partial}{\partial s^i}$ for $i = 1,\ldots,k$. If $S\subseteq W$ is an embedded codimension-$k$ submanifold and $p$ is a point of $s$ such that $T_pS$ is complementary to the span of $(V_1\vert_p,\ldots,V_k\vert_p)$, then the coordinates can also be chosen such that $S\cap U$ is the slice defined by $s^1=\cdots=s^k=0$.

At $p = (1,0,0)$ we have that $\text{span}\left\{X\vert_p,Y\vert_p\right\} = \text{span}\left\{\frac{\partial}{\partial y},\frac{\partial}{\partial z}\right\}$, so if we take $S$ to be the $x$-axis, then we can apply the theorem. However, since $S$ is not just a point, I am having trouble finding the explicit parameterization $(x(s,t), y(s,t), z(s,t))$. Agian, any help on how I need to go about this is greatly appreciated. Thank you!
 A: One can do this without directly thinking about or computing flows.
Hint Taking a general $1$-form $$\omega_0 := p(x, y, z) \,dx + q(x, y, z) \,dy + r(x, y, z) \,dz$$ on $\Bbb R^3$, imposing $\omega(X) = \omega(Y) = 0$, and solving for the coefficient functions $p, q, r$, shows that the annihilator of $X$ and $Y$ at each point is spanned by $$\omega_0 := dx - z \,dy - y \,dz,$$
but this form is exact:
$$\omega_0 = d(x - y z) .$$
(Involutivity of the distribution $D := \operatorname{span}\{X, Y\}$ implies that any annihilator $\omega$ of $D$ satisfies $d\omega = \theta \wedge \omega$ for some $1$-form $\theta$, and since $\Bbb R^3$ is simply connected, if $\omega$ vanishes nowhere there is a function $f$ such that $f \omega$ is exact.)

 So, by construction, $D$ is tangent everywhere to the level sets $G^{-1}(c)$, where $G(x, y, z) := x - y z$. In particular, $G(p) = G(1, 0, 0) = 1$, so the (maximal) integral submanifold $S$ of $D$ through $p$ is $S = G^{-1}(1) = \{x - y z = 1\}$. Rearranging the defining equation gives $x - y z - 1 = 0$, so we can take $$F(x, y, z) = x - y z - 1 .$$ Rearranging again to solve for $x$ gives $x = y z + 1$, giving the global parameterization $$(y, z) \mapsto (y z + 1, y, z)$$ of $S$.

A: With $g (x,y,z)=x-yz-1$, I think that it works. (To find it, I considered the equation that $Ker (g (x,y,z)) $ must satisfy, to find a candidate for $dg $ and then some $g $ that could work, and then check that it indeeds work). Considering that $[X,Y]=0$ over $\mathbb {R}^3$, we have  the unicity from the Frobenius theorem, but I'm not sure about that last point, I'm not really on point around that theorem.
Tell me if you think something there isn't right!
A: Ok, I finally figured out what I was doing wrong: I messed up the computations for the flows of the vector fields. The flows $\theta_t$ of $X$ and $\phi_t$ of $Y$ should be
\begin{align*}
\theta_t(x,y,z) &= (x + zt, y + t, z)\\
\phi_s(x,y,z) &= (x + ys, y, z + s).
\end{align*}
This would then give us that
$$\theta_t\circ \phi_s = \phi_s \circ \theta_t = (x + ys + zt + st, y + t, z + s).$$
If we take $S$ to then be the $x$-axis parameterized by $S(r) = (r, 0, 0)$, we can then use the given theorem to construct the function
$$\Psi_r(t,s) = \theta_t\circ \phi_s(r,0,0) = (r + st, t, s) = (x(s,t),y(s,t),z(s,t))$$
and it is fairly easy to check that
\begin{align*}
  \frac{\partial}{\partial t}\Psi_r(t,s) &= s\frac{\partial}{\partial x} + \frac{\partial}{\partial y} = z(s,t)\frac{\partial}{\partial x} + \frac{\partial}{\partial y} = X\\
  \frac{\partial}{\partial s}\Psi_r(t,s) &= t\frac{\partial}{\partial x} + \frac{\partial}{\partial z} = y(s,t)\frac{\partial}{\partial x} + \frac{\partial}{\partial z} = Y.
\end{align*}
Then if we want to construct the zero set $F(x,y,z) = 0$ that coincides with this surface passing through $(1,0,0)$, we simply pick $r = 1$ and solve the rest to get that $F(x,y,z) = x - yz - 1$.
