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}
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.