Solving hyperbolic system of pde using matlab As part of my MPhil research work am to solve the system of nonlinear PDEs below using matlab. It's my first time working with matlab and I am finding it difficult generating the code to solve the problem. I will be very glad if anyone can help me. The initial conditions are
\begin{equation}
(u_0,T_0,\rho_0,p_0)=(80m/s,290K,0.83kg/m^3,2001.23Pa)
\end{equation}
\begin{equation}
\begin{pmatrix}
\rho\\
\rho u\\
\epsilon \rho
\end{pmatrix}_t
+
\begin{pmatrix}
\rho u\\
\rho u^2+p\\
(\epsilon \rho +p)u
\end{pmatrix}_x
-
\begin{pmatrix}
0\\
\rho g \sin(\theta)\\
\rho gu \sin(\theta)
\end{pmatrix}
=\begin{pmatrix}
0\\
0\\
\frac{4h}{D}(3003-T)
\end{pmatrix}
\end{equation}
where 
$p=\frac{\rho RT}{1-\rho b}-a\rho^2$ and $\epsilon=c_vT-a\rho+\frac{u^2}{2}$
other constants are 
\begin{equation}
a=0.23, b=0.000043, h=100, D=0.5, g=10m/s^2,\theta=30^o, c_v=2.23, R=8.3145 
\end{equation}
 A: This is a 1D hyperbolic system of partial differential equations with relaxation
$$
\frac{\partial}{\partial t} \boldsymbol{U} + \frac{\partial}{\partial t} \boldsymbol{F}(\boldsymbol{U}) = \boldsymbol{R}(\boldsymbol{U}) \, ,
$$
with unknowns $\boldsymbol{U} = (\rho,\rho u,\epsilon \rho)^\top$, physical flux $\boldsymbol{F}(\boldsymbol{U}) = (\rho u,\rho u^2 + p,(\epsilon \rho + p)u)^\top$, and relaxation function $\boldsymbol{R}(\boldsymbol{U}) = (0, \rho g \sin\theta , \rho g u\sin\theta + \frac{4h}{D}(3003-T))^\top$. For this kind of problem, a convenient strategy is splitting, which results in a propagation part
$$
\frac{\partial}{\partial t} \boldsymbol{U} + \frac{\partial}{\partial t} \boldsymbol{F}(\boldsymbol{U}) = \boldsymbol{0}
$$
and a relaxation part
$$
\frac{\partial}{\partial t} \boldsymbol{U} = \boldsymbol{R}(\boldsymbol{U}) \, .
$$
Numerically, both parts can be solved successively at each time step -- e.g. perform Strang splitting -- with a dedicated numerical method. For instance, the propagation part is solved using the Lax-Friedrichs scheme, and the relaxation part is solved using the backward Euler method.
