I need to solve the following system of quasi-linear partial differential equations, where the unknowns $f(x,t), g(x,t)$ are smooth functions on $\mathbb R^2$,


$$ \dfrac{\partial}{\partial t}f=\dfrac{\partial}{\partial x}((f^2+g)f),$$
$$ \dfrac{\partial}{\partial t} g=\dfrac{\partial}{\partial x}((f^2+g)g),$$


I think I am supposed to use the method of characteristics to find a characteristic surface, but I don't know how to do this.


(I asked a more or less equivalent question on Math Stackexchange [here](https://math.stackexchange.com/questions/2397469/using-method-of-characteristics-to-solve-a-system-of-quasilinear-pdes), and got no answers)



Using method of characteristics to solve a system of first order quasilinear PDEs
The method of characteristics will tell you that, for any given solution $\bigl(f(x,t),g(x,t)\bigr)$, the characteristic curves are given by the foliations
$$
\mathrm{d}x + \bigl(3f(x,t)^2{+}2g(x,t)\bigr)\,\mathrm{d}t = 0
\quad\text{and}\quad
\mathrm{d}x + \bigl(f(x,t)^2{+}g(x,t)\bigr)\,\mathrm{d}t = 0.
$$
Thus, the general solutions of the system are hyperbolic, except when the characteristics are confounded (i.e., double), which is when $2f(x,t)^2+g(x,t) = 0$.


One way to see this is to note that the above equations are equivalent to the *conservation laws*
$$
\mathrm{d}\bigl(f\,\mathrm{d}x + (f^2+g)f\,\mathrm{d}t\bigr) = 
\mathrm{d}\bigl(g\,\mathrm{d}x + (f^2+g)g\,\mathrm{d}t\bigr) = 0,
$$
which is equivalent (away from the double characteristic case) to the following equations expressed as decomposable $2$-forms:
$$
d(f^2{+}g) \wedge(\mathrm{d}x + (3f^2{+}2g)\,\mathrm{d}t) = 
(g\,\mathrm{d}f - f\,\mathrm{d}g) \wedge \bigl(\mathrm{d}x + (f^2{+}g)\,\mathrm{d}t\bigr) = 0
$$ 


Now, the only solutions that have double characteristics everywhere belong to the $1$-parameter family $(f,g) = (c,-2c^2)$ for some constant $c$. You can see this if you substitute $g = -2f^2$ in the above equations, which then imply that $f\_t=0$ and $(f^3)\_x=0$, so $f$ (and hence $g$) must be constant.


The cases of simple characteristics satisfying an integrable relation are when either $\mathrm{d}(f^2+g) = 0$, which is the case that Robert Israel solved, or when $g\,\mathrm{d}f - f\,\mathrm{d}g = 0$,
which is the case that $f = c\_1 u(x,t)$ and $g = c\_2 u(x,t)$ for some function $u$ where $c\_1$ and $c\_2$ are constants. In this case, the function $u(x,t)$ will be given implicitly as a function of $x$ and $t$ by an equation
$$
x + t\,(3{c\_1}^2u^2+2c\_2u) - U(u) = 0
$$
where $U$ is an arbitrary function of $u$ subject to the condition that the equation above defines $u$ as a smooth function of $x$ and $t$. [NB: It's actually slightly more general than this. What is true is that, for a solution, the functions $u$ and $x + t\,(3{c\_1}^2u^2+2c\_2u)$ have the same level sets, i.e., they are functionally related. However, assuming that $x + t\,(3{c\_1}^2u^2+2c\_2u)$ is a function of $u$ covers most cases of interest.]


For the general solution, you'll need to study the characteristic distributions determined by the initial conditions. What you can say is that, if you want to specify initial conditions along a curve $x=x\_0(s)$, $t=t\_0(s)$ in the $xt$-plane in the form 
$$
f\bigl(x\_0(s),t\_0(s)\bigr) = f\_0(s)\quad\text{and}\quad
g\bigl(x\_0(s),t\_0(s)\bigr) = g\_0(s),
$$
then there will be a unique smooth solution in a neighborhood of this curve as long as the initial data is both *noncharacteristic*, i.e., as long as 
$$
x'\_0(s) + \bigl(3f\_0(s)^2+2g\_0(s)\bigr) t'\_0(s) \not=0\quad\text{and}\quad
x'\_0(s) + \bigl(f\_0(s)^2+g\_0(s)\bigr) t'\_0(s) \not=0,
$$
and *hyperbolic*, i.e., as long as $g\_0(s)+2f\_0(s)^2\not=0$. Typically, the initial conditions will only allow a smooth solution in a neighborhood of the initial curve, with shocks developing at finite distance from the curve. In that case, if you want to continue the solution past the shocks, you'll need to study the Rankine-Hugoniot jump conditions. It might be helpful to look at the classic work of James Glimm on the existence of hyperbolic solutions in the quasi-linear case with conservation laws.


It's possible that the full system is integrable by the method of Darboux (generalizing the method of Riemann invariants), but I haven't checked for that. 


