I'm preparing for an exam by looking at an old exam archive. I encountered this problem and am not sure how to proceed after showing that the distribution is integrable. I haven't learned about differential forms yet. Any hints would be wonderful. The problem is:
"Consider the distribution in $\mathbb{R}^3$ spanned by the two vector fields
$$V = \partial_x + 2xy \partial_z, \hspace{10mm} W = x \partial_x + \partial_y + (2x^2y + x^2 - 2y) \partial_z.$$
Show that this distribution is integrable and find an explicit formula for the integral submanifold passing through the point $(0,0,z_0)$."
To see that this distribution, call it $\Delta$, is integrable, we want to show that $[\Delta, \Delta] \subseteq \Delta$. We may check the Lie bracket on just the generating vector fields. In this case, there is only the one $[V,W] = \partial_x + (4xy + 2x)\partial_z-2xy\partial_z-2x\partial_z=\partial_x + 2xy\partial_z = V$. So it is integrable by Frobenius' Theorem.
Now, we want to find the integrable submanifold. I think this involves solving a PDE in the end and perhaps choosing nice coordinates so that the formula has the form $(x,y,f(x,y))$. I'm very weak in PDE's however; well, I haven't learned them at all. One idea is that perhaps we could consider $V' = V$ and $W' = W-xV=\partial_y +(x^2-2y)\partial_z$ instead. It seems slightly easier to work with.
Differential forms are the way to go. You can check that this distribution is given by $\omega = 0$ for $\omega = dz - 2xy\,dx + (2y-x^2)\,dy$. Note that $d\omega = 0$ (not just $0$ mod $\omega$). Indeed, you can see that $\omega = df$ for $f(x,y,z) = z+y^2-x^2y$. You want the level surface of $f$ passing through $(0,0,z_0)$, which is the surface $z+y^2-x^2y = z_0$.