So with some guidance from @TedShifrin I think I got it.
First since $\omega$ is decomposable and closed taking the exterior derivative, $0=d\phi \wedge \theta -(1) \phi \wedge d\theta$. Hence, $d\phi \wedge \theta = \phi \wedge d\theta$. Wedging the expression with $\phi$ gives $d\phi \wedge \theta \wedge \phi=0$, since the right side has $\phi \wedge \phi$ in there. Similarly, $d\theta \wedge \theta \wedge \phi =0 $.
Hence, $d\theta$ and $d\phi$ are linearly dependent on $\theta \wedge \phi$, hence $d\theta, d\phi \in \langle \theta, \phi \rangle $.
The above sentence is incorrect (I have not figured out how to cross it out. I would like to leave it there for instructional purposes).
Rather, I need to prove the following general identity:
Let $\omega^1, ... \omega^k$ be locally defining independent one forms for a $n-k$ rank distribution $D$. Then $D$ is involutive iff $d\omega^i \wedge \omega^1 \wedge ... \wedge \omega^k =0$ for all $1\leq i \leq k$. One direction is standard by the $1$-form criterion for involutivity, which says that $D$ is involutive iff for any 1-form $\eta$ annihilating $D$, $d\eta$ annihilates $D$.
Conversely, suppose the above holds. Take the independent $\omega^1, ...\omega^k$ and extend to a smooth local coframe $\{ \omega^1, ...\omega^k, \alpha^{k+1}, ... \alpha^n \}$ with a dual frame $\{ E_1, ... E_k, A_{k+1},...A_n \}$. Note the $A_j$'s span $D$ Then $d\omega^{i}=\Sigma b_{ij}\omega^{i}\wedge \omega^{j} + \Sigma c_{kl} \omega^{k} \wedge \alpha^{l} + \Sigma d_{xy}\alpha^{x}\wedge \alpha^{y}$ with sums over appropriate increasing indices.
Now, wedge this with $\omega^1 \wedge ... \wedge \omega^k$ and evaluate at $(A_x, A_y, E_1, ..., E_k)$. This then shows $d_{xy}=0$. Hence in fact $d\omega^{i}=\Sigma b_{ij}\omega^{i}\wedge \omega^{j} + \Sigma c_{kl} \omega^{k} \wedge \alpha^{l}$ which when evaluated at any $(A_x, A_y)$ is equal to zero. Hence by the 1-form criterion $D$ is involutive.
Now the next paragraph can in fact be skipped.
Thus, if $\lambda=\alpha \wedge \theta + \beta \wedge \phi$, we have $d\lambda=d\alpha \wedge \theta - \alpha \wedge d\theta + d\beta \wedge \phi - \beta \wedge d\phi \in \langle \theta, \phi \rangle$, hence the ideal is a differential ideal.
Thus $D$ defines and involutive distribution, which is integrable by the Frobenius theorem.
Thus there is a flat chart where the distribution, call it $D$ is spanned by $\frac{\partial}{\partial x^{3}}, ... , \frac{\partial}{\partial x^{n}}$ and annihilated by the defining $1$-forms $\theta, \phi$. The only way for this to happen is for $\theta = f_{1}dx^1 + f_{2}dx^2$ and $\phi = g_{1}dx^1 + g_{2}dx^2$ Taking the wedge, we get $\omega =(g_{1}f_{2}-g_{2}f_{1})dx^1 \wedge dx^2=F(x^1, ..., x^n)dx^1\wedge dx^2$.
Now let us take the exterior derivative of $\omega.$
$d\omega=\Sigma \frac{\partial F}{\partial x^k}dx^k \wedge dx^1 \wedge dx^2=0$ implying the partial derivatives vanish. Hence $F=F(x^1, x^2)$
Now define a new coordinate system by $y^1=\int_{0}^{x^1}F(t, x^2)dt, y^2=x^2, ... , y^n=x^n$. The Jacobian of this transform is non-zero. Further $dy^1=Fdx^1+ \frac{\partial}{\partial x^2}(y^1=\int_{0}^{x^1}F(t, x^2)dt)dx^2.$ Hence in this coordinate system $\omega=Fdx^1 \wedge dx^2=dy^1 \wedge dy^2$.