Note that there is an easy integrating factor to guess for your $\omega$: Multiply through by $f=\dfrac1{(x^2+y^2)z}$ and the resulting form is exact. Indeed, $f\omega = d(-\log r + \log z)$, where $r=\sqrt{x^2+y^2}$. Integral manifolds are then $z=cr$ for various constants $c$.
Following Sophus Lie, there is a good way to look for such integrating factors. If you see an obvious one-parameter group leaving the distribution invariant (or sending $\omega$ to a functional multiple of $\omega$), let $V$ be its infinitesimal generator and let $f=1/\omega(V)$. Then you can check that $f\omega$ will be closed. In your case, the obvious one-parameter group is given by radial stretches, and $V=x\dfrac{\partial}{\partial x}+y\dfrac{\partial}{\partial y}+z\dfrac{\partial}{\partial z}$. Unfortunately, the numbers work out to $\omega(V)=0$ — in other words, $\mathscr L_V\omega = 0$ rather than just a general functional multiple of $\omega$. But, still, we will see the appearance of $(x^2+y^2)z$ when we do this computation.
EDIT: A natural way to avoid the issue with $z=0$ about which you're worried is to introduce a homogenizing variable $u$ with $z=ur$, $r=\sqrt{x^2+y^2}$. Then
$$\omega = -zr\,dr + r^2\,dz = -ur^2\,dr + r^2(u\,dr + r\,du) = r^3\,du,$$
and so, naturally, $\dfrac 1{r^3}\omega = du$ is exact. (So, here, in the new coordinate system, $1/r^3$ does become the integrating factor!) Oh, and finally, the integral manifolds are patently obvious: They're given by $u=c$ for all constants $c$.
Expanding $(\lambda d\alpha + d\lambda \wedge \alpha)^n$, you should get
$$
\sum_{k=0}^n c_k\cdot (\lambda d\alpha)^{n-k} \wedge (d\lambda \wedge \alpha)^k
$$
where $c_k$ is a combinatorial constant.
Notice that $c_0=1$.
Now, any term with $k>0$ has some $\alpha$ term inside, and after a suitable permutation in the wedge products, one can write
$$
(\lambda d\alpha + d\lambda \wedge \alpha)^n = \lambda^ n (d\alpha)^n + \alpha \wedge \beta,
$$
for $\beta$ a $(2n-1)$-form.
Finally, $\alpha$ being a $1$-form, we have $\alpha\wedge \alpha = 0$, and therefore
$$
\lambda \alpha \wedge (d(\lambda \alpha))^n = \lambda^{n+1}\alpha \wedge (d\alpha)^n + \lambda \underbrace{\alpha \wedge \alpha}_{=0} \wedge \beta = \lambda^{n+1} \alpha\wedge (d\alpha)^n.
$$
Best Answer
Let's think for a moment about what it means to say $\omega\wedge d\omega = 0$ when $\omega\in\Omega^1(M)$. If you consider extending $\omega$ to a basis $\{\omega=\omega_1,\omega_2,\dots,\omega_n\}$ for the $1$-forms locally, you can easily convince yourself that $d\omega = \omega\wedge\eta$ for some $1$-form $\eta$, again locally. Then (locally) it is clear that $d\omega(X,Y) = 0$ whenever $\omega(X)=\omega(Y)=0$. Thus, $\omega\wedge d\omega = 0$ implies involutivity.
On the other hand, if the distribution is involutive, then you already know that $d\omega(X,Y) = 0$ whenever $\omega(X)=\omega(Y)=0$. But then this tells us that $d\omega = \omega\wedge\eta$ for some $1$-form $\eta$; for, if not, we would have $$d\omega=\sum\limits_{i=2\\i<j}^n f_{ij}\omega_i\wedge\omega_j$$ for some functions $f_{ij}$, with some coefficient $f_{i'j'}$ nonzero (at a particular point). Choosing $X,Y$ so that $\omega_{i'}(X)=1$ and $\omega_k(X)=0$ for all $k\ne i'$ (including $k=1$) and $\omega_{j'}(Y)=1$ and $\omega_k(Y)=0$ for all $k\ne j'$ (likewise) shows that $d\omega(X,Y)\ne 0$.