Let $a,b$ two smooth functions from the open square $I^{2}$ in $\mathbb{R}^{2}$ to $\mathbb{R}^{4}$. In particular, assume $a(t,u)$ and $b(t,u)$ be linearly independent for all $(t,u) \in I^{2}$.


I need to study the PDE problem in $ x \in C^{\infty}(I^{2},\mathbb{R}^{4})$ given by the underdetermined system
$$
\begin{cases}
\displaystyle \frac{\partial x}{\partial t} \cdot a = 0 \\
\displaystyle \frac{\partial x}{\partial u} \cdot b = 0 \\
\displaystyle \frac{\partial x}{\partial t} \cdot b = \frac{\partial x}{\partial u} \cdot a
\end{cases}
$$
and the auxiliary condition
$$x(\cdot,0)= \frac{\partial x}{\partial u}(\cdot,0)=0\,.$$


**Question**: how would you go about solving this problem, i.e. representing its solution set?


Before somebody legitimately asks for motivation, let me finally say that this problem is linked to the isometric flexibility of a surface in $\mathbb{R}^{4}$.



Underdetermined system of linear PDEs
Is there anything else that you are not telling us about $a$ and $b$? The particulars of these two vector-valued functions have a great influence on what the general solution of the system 
$$
a\cdot x\_t = b\cdot x\_u = a\cdot x\_u - b\cdot x\_t = 0\tag 1
$$
looks like.


For example, take the very special case in which $a$ and $b$ are *constant*. Then the above equations become 
$$
(a\cdot x)\_t = (b\cdot x)\_u = (a\cdot x)\_u - (b\cdot x)\_t = 0\tag 2
$$


Let $a^\*$ and $b^\*$ be the two linear combinations of $a$ and $b$ that satisfy
$$
a\cdot a^\* = b\cdot b^\* = 1\qquad\text{while}\qquad a\cdot b^\* = b\cdot a^\* = 0,
$$
and write $x = v\,a^\* + w\,b^\* + c$ where $c:I^2\to\mathbb{R}^4$ satisfies $a\cdot c = b\cdot c = 0$. Then the above equations become $v\_t = w\_u = v\_u-w\_t=0$, which implies that there are constants $v\_0$, $w\_0$, and $r\_0$ such that
$$
x = (v\_0 + r\_0\, u)\,a^\* + (w\_0 + r\_0\, t)\,b^\* + c.\tag 3
$$
Moreover, $c$ is arbitrary subject to the equations $a\cdot c = b\cdot c = 0$. This is the general solution of the system in the case that $a$ and $b$ are (linearly independent) constants. Imposing the given 'auxilliary conditions' forces $v\_0=r\_0=w\_0=0$ and $c(t,0)=c\_u(t,0)=0$, so the solution reduces to $x(t,u) = u^2w(t,u)$ where $w:I^2\to\mathbb{R}^4$ is any mapping satisfying the two linear constraints $a\cdot w = b\cdot w = 0$. Thus, the general solution subject to the auxilliary conditions still depends on two arbitrary functions of two variables. (A similar story holds whenever $a$ and $b$ take values in a fixed $2$-dimensional subspace $S\subset\mathbb{R}^4$.)


At the other extreme, consider the case in which the $5$ vector-valued functions $a$, $b$, $a\_t$, $b\_u$, and $a\_u{-}b\_t$ span $\mathbb{R}^4$ at every $(t,u)\in I^2$. Let us write
$$
p\,a\_t + q\,b\_u + r\,(a\_u{-}b\_t) + f\,a + g\,b = 0\tag 4
$$
for some functions $p$, $q$, $r$, $f$ and $g$, where $p$, $q$, and $r$ do not simultaneously vanish. Up to a common nonzero multiple, these $5$ functions are determined by $a$ and $b$, so they can be regarded as known.


Set 
$$x\cdot a = v\qquad \text{and} \qquad x\cdot b = w.\tag 5
$$  

Then, with the given equations, we have the three identities
$$
x\cdot a\_t = v\_t\qquad x\cdot b\_u = w\_u\qquad x\cdot(a\_u-b\_t) = v\_u-w\_t\,.\tag 6 
$$
Given the relation (4), we then get a single linear first order PDE for the pair $(v,w)$:
$$
p\,v\_t + q\,w\_u + r\,(v\_u-w\_t) + f v + g\,w = 0.\tag 7
$$
Given a solution $(v,w)$ to (7), the five inhomogeneous linear equations (5)$+$(6) for $x$ are consistent and uniquely determine a solution $x$ to the original system of $3$ equations. Thus, the solutions of the original overdetermined system of three equations for four unknowns have now been expressed in terms of the solutions to a single linear PDE for two unknowns. Imposing the 'auxilliary conditions' is equivalent to requiring that the five functions $v$, $w$, $v\_t$, $w\_u$, and $v\_u{-}w\_t$ and their first derivatives with respect to $u$ should all vanish along the edge $u=0$ in $I^2$.


Methods for writing down the solutions of an equation such as (7) can be found in the 19th century PDE literature. I can say more about this if you are interested. The most elementary thing to do is to simply take one of the pair $(v,w)$ to be an arbitrary function on $I^2$ and then use the integrating factor method to solve the resulting linear PDE for the other. 


There are, of course, more sophisticated things one can do to reduce the problem further to an 'explicit solution' without any integration involved. In particular, it can be shown that, for 'generic' $a$ and $b$ (i.e., for $(a,b)$ in a dense, open set in the $C^\infty$-topology on pairs of smooth maps $I^2\to\mathbb{R}^4$), there exist second-order linear differential operators $A$ and $B$ on $C^\infty(I^2)$ such that $(v,w) = \bigl(A(z),B(z)\bigr)$ solves (7) for every smooth function $z$ on $I^2$ and, moreover, every smooth solution $(v,w)$ of (7) is of the form $(v,w) = \bigl(A(z),B(z)\bigr)$ for some smooth function $z$ on $I^2$. The auxilliary conditions then translate into a pair of linear conditions on $z$ and its partials up to fourth order along the edge $u=0$.


In between these two extreme cases of hypotheses on $a$ and $b$, there are some exceptional cases that require separate handling. Also, I can say more about this, if you want.


