Method of characteristics for systems of PDE (vs. Lewy's example) Main question: How does the method of characteristics generalize for systems of first order PDE, as opposed to scalar PDE?  Namely, is there such a generalization at all, and if so what information can it provide (e.g. implicit solution formulas, local existence, etc.)?  Also, are there any references for this matter?
Corollary question/motivation: The method of characteristics provides a local existence result for scalar, first order PDE.  Lewy's example of a PDE with no solutions on any open set (as presented on wiki example) is also first order.  Presumably the failure of the aforementioned local existence result here is due to the fact that Lewy's example is effectively a system consisting of the real and imaginary parts of u, correct?  So any generalization of the method of characteristics to systems must be rather less effective than the scalar version...
Further motivation: The method of characteristics provides a very useful way to solve scalar conservation laws and Hamilton-Jacobi equations numerically.  If there is a generalization of the method of characteristics to systems, does it provide a similarly powerful tools for numerically solving systems of PDE?
 A: Any partial differential equation can be written as a first order system, so basically what you are asking is how to generalize method of characteristics to general PDEs and systems. As you have already observed, in general it is hopeless. However, there are some things that can be said.
For general PDEs and systems, the notion of characteristic surfaces plays a crucial role, which can be considered as a substitute for characteristic curves. Further, when we study high frequency asymptotics of (or how singularities propagate under) a general linear PDE, we are led to a fully nonlinear first order equation (of Hamilton-Jacobi type), which can be solved by the method of characteristics. The "characteristic curves" that arise are called bicharacteristics, which of course lie in a higher dimensional space (with double the dimension of the space of independent variables). Bicharacteristics refine the notion of characteristic surfaces, and are especially important for hyperbolic and dispersive equations, and there are many numerical methods based on this structure. Basically, what this means is that to solve the wave equation, you start with geometric optics, or to solve the Schrödinger equation, you start with classical mechanics. On the analytic side this idea eventually led to the powerful technique of Fourier integral operators.
Another thing I wanted to mention is that in 1 space dimension, if your (hyperbolic) system is simple enough, locally it would look like a collection of first order scalar equations, and you can use method of characteristics componentwise.
