Solution (source code)

= Solution

Write the system as $U_x+A(U)U_y=0$ with $U=(u,v)^T$ and
$$
A=\begin{pmatrix}v&u\\1&v\end{pmatrix},\qquad\det(A-\lambda I)=(v-\lambda)^2-u.
$$
Its two <characteristic speeds> are $\lambda_\pm=v\pm\sqrt u$, so it is \b[strictly hyperbolic precisely when $u>0$]. At $u=0$ they coincide and the matrix is not diagonalizable; for $u<0$ the speeds are nonreal.

For $u>0$ define $R_\pm=v\pm2\sqrt u$ and $D_\pm=\partial_x+(v\pm\sqrt u)\partial_y$. Substitution from the two differential equations gives
$$
\begin{aligned}
D_\pm R_\pm
&=v_x+(v\pm\sqrt u)v_y
\ \pm u^{-1/2}[u_x+(v\pm\sqrt u)u_y]\\
&=(-u_y\pm\sqrt u\,v_y)\pm u^{-1/2}(\pm\sqrt u\,u_y-uv_y)=0.
\end{aligned}
$$
Therefore
$$
\boxed{v\pm2\sqrt u\text{ is constant on curves }\frac{dy}{dx}=v\pm\sqrt u.}
$$
These are the two <Riemann invariants>; their constancy applies within smooth regions, not across a <shock>.