Explain why the characteristic curves generate the solution surface I recently came across this problem in one of my differential equation texts:

Explain why the characteristic curves generate the solution surface of $P(x,y) \frac{\partial z}{\partial x} + Q(x,y) \frac{\partial z}{\partial y} = R(x,y,z)$ assuming $P,Q$ are Lipschitz in $x,y$ and $R$ is Lipschitz in $z$.

Now I know the proof of this fact, but I think this is not the desired answer. Is there perhaps a (geometrical) explanation of why this might be the case?
I think it might be possible to reason as follows: One can show that $(P,Q,R)$ lies in the tangent plain of any point on the solution surface. $(P,Q,R)$ is by definition also tangent to any point on the characteristic curves...
 A: Let's look at the proof of the fact and try to extract the geometric meaning in its steps. We are trying to solve:
$$
P(x, y)u_x + Q(x, y) u_y - R(x, y, u) = 0 \label{1}\tag{1}.
$$
The characteristic equations are
$$
\begin{cases}
(a) \ \ (\dot x(s), \dot y (s)) = (P(x, y), Q(x, y)) \\
(b) \ \ \dot z(s) = R(x, y, z).
\end{cases}
$$
We can write this as just a single equation for the triple $(x(s), y(s), z(s)) = \gamma(s)$:
$$
\dot \gamma(s) = (P(x, y), Q(x, y), R(x, y, z))
$$
(note that all arguments on the right hand side depend on the parameter $s$). These are our characteristic curves.
Now we can argue as you say. The point is exactly the tangency you mentioned. The solution surface is the graph of a solution $z = u(x, y)$ in $(x, y, z)$-space. The first step is to recast the PDE \eqref{1} into a geometric interpretation. For this we use that a normal vector to a graph of a function $z = f(x, y)$ is given by $(f_x, f_y, -1)$. So the normal vector to the graph of $z = u(x, y)$ is
$$
(u_x, u_y, -1).
$$
If we dot product this with $\vec{V} = (P, Q, R)$, as you say, we see that
$$
\vec{V}\cdot (u_x, u_y, -1) = Pu_x + Qu_y - R = 0
$$
by the PDE \eqref{1}. So, saying that $u$ solves \eqref{1} is equivalent to saying that $\vec{V}$ is tangent to the graph of $u$. Now let's look at the characteristic curves. The next point is that the characteristic curves are integral curves of $\vec{V}$. This we have actually already shown above: $\gamma$ satisfies
$$
\dot \gamma(s) = (P, Q, R) = \vec{V}(\gamma(s)).
$$
The last equality just follows by rewriting things in terms of our new $\vec{V}$ notation.
Thus, the curves $\gamma$ must lie entirely within the solution surface, since they are integral curves of the vector field $\vec{V}$ which, we computed above, is tangent to the solution surface. As for the question of why they must cover the whole solution surface: given any $(x_0, y_0, z_0 = u(x_0, y_0))$ on the solution surface, we have the vector field $\vec{V}$ here and can solve the characteristic ODE $\dot\gamma = \vec{V}(\gamma)$ with the initial condition $\gamma(0) = (x_0, y_0, z_0)$.
A: Part (1) : Given a solution to the PDE, denoted by $u(x,y)$, and the surface generated by its graph in ${\Bbb R}^3$ it is clear that a trajectory of the corresponding ODE: $$\dot{x_t}=P(x_t,y_t), \dot{y_t}=Q(x_t,y_t), \dot{z_t}=R(x_t,y_t,z_t),$$ starting from a point in the surface will stay in the surface. This follows from
$$ \frac{d}{dt} u(x_t,y_t) = \dot{x_t} \frac{\partial u}{\partial x} + \dot{y_t} \frac{\partial u}{\partial y} = P(x_t,y_t)\frac{\partial u}{\partial x} + Q(x_t,y_t) \frac{\partial u}{\partial y} =R(x_t,y_t,u(x_t,y_t)).
$$
As $\dot{z_t} = R(x_t,y_t,z_t)$ we see that $t\mapsto z_t$ and $t\mapsto u(x_t,y_t)$ verify the same ode, so starting from the same point, they
must be identical. Now, this is the easy direction. You may interpret in terms of tangent vectors being orthogonal to surface-normals etc, but this is perhaps not so relevant.
A much harder part, and perhaps the intention with the posed question, is to see why a collection of solution curves generates a surface satisfying the PDE and having a natural geometric interpretation. So let me attack this part. We assume that $P,Q,R$ are at least $C^1$:
Part (2): The vector field $V=(P(x,y),Q(x,y))$ in the plane gives rise to a local flow $\phi^t$ in ${\Bbb R}^2$, given by solving $\dot{x}=P$, $\dot{y}=Q$ with initial conditions $(x(0),y(0))= \xi_0=(x_0,y_0)$.
If, in addition we are given some initial value $z(0)=z_0$, we may
integrate the ode $\dot{z}=  R(x(t),y(t),z(t))$, obtaining a curve $(x(t),y(t),z(t))$ in 3-space. We say that this curve is a lift of the trajectory $(x(t),y(t))$ in the plane.
At first there seems to be no relation to the PDE, since
we have only produced some function $z(t)$ together with the flow-line in the plane. In order to get a (locally defined) surface we must assume that $V$ is non-zero at some given point $\xi_0$.  We then pick an arbitrary locally defined curve $\gamma$, transverse to the vector field $V$ at $\xi_0$. Let us write $\xi_0(s)=(x_0(s),y_0(s))$, $s\in I$ (with $I$ some small interval) for a $C^1$-regular parametrization of $\gamma$. When acting with the flow on $\gamma$ over some (small) time-interval $J$ we will sweep out a domain  in the $(x,y)$-plane. We write it as: $$A = \{(x(s,t),y(s,t)) = \phi^t(\xi_0(s)) : s\in I, t\in J\}.$$
By the parameter dependency part of Cauchy-Lipshitz's theorem,
everything is at least $C^1$.
And since we required transversality, the map $(s,t)\mapsto (x(s,t),y(s,t))$ must be a local diffeo, meaning that we may use either $(s,t)$ or $(x,y)$ as local coordinates.
Now, pick in addition any $C^1$ function $z_0(s)$ along this curve and transport this curve along the flow using the ODE for $z$ discussed above.
The curve in ${\Bbb R}^3$ this time sweeps out a surface
$$ S = \{(x(s,t),y(s,t),z(s,t)) : (s,t)\in I\times J\},$$
i.e. a graph over the domain $A$. As just mentioned we may use local coordinates $(x,y)$ and identify  $u(x,y) \equiv {z}(s,t)$ at equivalent points. Again everything is $C^1$. The (Lie-) derivative of $u$ in the direction of $V$ and evaluated at a given point may then be calculated in two different ways by deftly joggling between the two local coordinates:
$$ L_V u  = P (x,y)\frac{\partial u}{\partial x} + Q (x,y)\frac{\partial u}{\partial y} = \frac{\partial x(s,t)}{\partial t} \frac{\partial u}{\partial x} + \frac{\partial y(s,t)}{\partial t} \frac{\partial u}{\partial y} = \frac{\partial z(s,t)}{\partial t} = R(x,y,u(x,y)).$$
By construction, this yields a $C^1$ solution to our PDE. Furthermore, as is readily seen  any (local) $C^1$-solution may be described in this way at points where $V$ is non-vanishing. At points where $V$ vanishes the method of characteristics does not apply.
