I have the vector fields $v_{1} = x \partial_y - y \partial_x + z \partial_w - w \partial_z$ and $v_{2} = z \partial_x - x \partial_z + w \partial_y - y \partial_w$ on $S^{3} \subset \mathbb{R}^4$.
I have shown that the Lie bracket is zero: $[v_1, v_2] = 0$ and by Frobenius' {$v_1$, $v_2$} form an integral system.
Now I am asked to find the integral submanifolds in $S^3$, but am unsure about how to do this.
I am currently working on the same problem. I believe the idea is to use something akin to the method of characteristics for partial differential equations.
Suppose that $$\gamma(t) = \begin{pmatrix}x(t)\\y(t)\\z(t)\\w(t)\end{pmatrix}$$ is a curve on $S^3$ whose tangent vectors correspond to the vector field $v_1$. Also suppose that $\gamma$ is parameterized so that $x(0) = x_0$, $y(0) = y_0$, $z(0) = z_0$, and $w(0) = w_0$. If this is the case, then the coordinate functions satisfy the following system of ordinary differential equations:
\begin{align} &\dfrac{dx}{dt} = -y\\[0.3cm] &\dfrac{dy}{dt} = x\\[0.3cm] &\dfrac{dz}{dt} = -w\\[0.3cm] &\dfrac{dw}{dt} = z\\ \\ & \text{with $x(0)=x_0$, $y(0) = y_0$, $z(0) = z_0$, and $w(0) = w_0$} \end{align}
Next notice the following:
\begin{align} & \dfrac{\dfrac{dx}{dt}}{\dfrac{dy}{dt}} = \dfrac{dx}{dy} = -\dfrac{y}{x }\\ \\ & \dfrac{\dfrac{dz}{dt}}{\dfrac{dw}{dt}} = \dfrac{dz}{dw} = -\dfrac{w}{z } \end{align}
Solving these ODE's, we arrive at the result: $$\gamma(t) = \begin{pmatrix}\sqrt{A} \cos t\\ \sqrt{A}\sin t\\ \sqrt{B}\cos t\\ \sqrt{B}\sin t\end{pmatrix}$$ Where $A$ and $B$ are determined from the initial conditions to be $x_0^2+y_0^2$ and $z_0^2+y_0^2$ respectively.
Next consider use the same method to flow out from the point $p_t=(x_t, y_t, z_t,w_t)$ by a different curve $\eta(s)$ that satisfies the initial condition $\eta(0) = p_t$ and whose tangent vectors correspond to $v_2$. Following this method should lead you to the complete solution.