To solve 2D quasilinear systems of conservation laws
$$
{\bf u}_t + {\bf A}({\bf u})\, {\bf u}_x + {\bf B}({\bf u})\, {\bf u}_y = {\bf 0}
$$
numerically,
various strategies can be followed:
- Implement a 2D finite-volume scheme, such as the 2D Lax-Friedrichs method
\begin{aligned}
{\bf u}_{i,j}^{n+1} &= \frac{{\bf u}_{i-1,j}^{n} + {\bf u}_{i+1,j}^{n} + {\bf u}_{i,j-1}^{n} + {\bf u}_{i,j+1}^{n}}{4} \\ & - {\bf A}({\bf u}_{i,j}^n) \frac{{\bf u}_{i+1,j}^{n} - {\bf u}_{i-1,j}^{n}}{2\, \Delta x/\Delta t}
-{\bf B}({\bf u}_{i,j}^n) \frac{{\bf u}_{i,j+1}^{n} - {\bf u}_{i,j-1}^{n}}{2\, \Delta y/\Delta t} .
\end{aligned}
- Use the Method of Lines, which consists in integrating the semi-discrete system
$$
\dot{\bf u}_{i,j} = -{\bf A}({\bf u}_{i,j})\, ({\bf u}_x)_{i,j} - {\bf B}({\bf u}_{i,j})\, ({\bf u}_y)_{i,j}
$$
in time.
- Introduce the dimensional splitting
\begin{aligned}
&{\bf u}_t + {\bf A}({\bf u})\, {\bf u}_x = {\bf 0} \\
&{\bf u}_t + {\bf B}({\bf u})\, {\bf u}_y = {\bf 0}
\end{aligned}
and integrate successively each part with 1D finite-volume schemes by using an operator splitting procedure.
For the one-dimensional scalar case, you will find several examples of MATLAB code on this site (e.g., (a), (b), (c)). Besides the representation of data in multi-dimensional arrays, the implementation of 2D methods is very similar.
(1) R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge university press, 2002. doi:10.1017/CBO9780511791253
(2) J.A. Trangenstein, Numerical Solution of Hyperbolic Partial Differential Equations, Cambridge university press, 2009. isbn:9780521877275
