Thanks for the clarification; I wasn't familiar with this terminology. I assume that the coefficients $a_i$, $b_i$, $c_i$, $d_i$, $e_i$, and $f_i$ are specified functions of $x_1,x_2,h_1,h_2$. (Let me know if this is not correct.)

Geometrically, what you have is a pair of $2$-forms on a $4$-dimensional manifold, which you have specified in coordinates $(x_1,x_2,h_1,h_2)$ as
$$
\begin{align}
\Upsilon_1 &= a_1\ dx_1\wedge dx_2 + b_1\ dh_1\wedge dx_2 + \cdots + e_1\ dh_2\wedge dx_2 + f_1\ dh_1\wedge dh_2\\\\
\Upsilon_2 &= a_2\ dx_1\wedge dx_2 + b_2\ dh_1\wedge dx_2 + \cdots + e_2\ dh_2\wedge dx_2 + f_2\ dh_1\wedge dh_2\ ,\\\\
\end{align}
$$
and you are looking for surfaces that are Lagrangian with respect to both $2$-forms, as these correspond to (generalized) solutions to your PDE system.

In order for this system to be 'symplectic', what you want is to find (or, better, determine whether there exist) combinations $\Omega = \lambda^1\ \Upsilon_1 + \lambda^2\ \Upsilon_2$ that satisfy $\Omega^2\not=0$ and $d\Omega =0$. Typically, this is 4 first order linear equations for the two unknown functions $\lambda^i$, so it's overdetermined, and there are 'curvature invariants' that determine when there are solutions and how many there are.

These are what are often called 'conservation laws' in some literature on this subject, and one place you could look (at least in the hyperbolic case; the elliptic case and parabolic cases are somewhat similar, though, at least in the formal setting) is in our paper with Phillip Griffiths and Lucas Hsu, *Hyperbolic exterior differential systems and their conservation laws, I & II*, which appeared in Selecta Mathematica (New Series) **1** (1995), pp. 21–112 and 265–323. (.dvi files of these papers are available at http://fds.duke.edu/db/aas/math/faculty/bryant/publications.html)