I have a system of 4 nonlinear 1st-order PDEs. I want to solve them numerically by method of lines, first discretizing space. This leads to the system of $N\times 4$ coupled ODEs.
$$ \mathbf{u}_{i} = (u^{1},u^{2},u^{3},u^{4})_{i},\quad i=1,\dots,N\\ \frac{\mathrm{d}u_{i}^{1}}{\mathrm{d}t} = q_{1}(\mathbf{u}_{i}) \\ \frac{\mathrm{d}u_{i}^{2}}{\mathrm{d}t} = q_{2}(\mathbf{u}_{i}) \\ \frac{\mathrm{d}u_{i}^{3}}{\mathrm{d}t} = F_{1}(\mathbf{u}_{i-k},\dots,\mathbf{u}_{i+l}) \\ \frac{\mathrm{d}u_{i}^{4}}{\mathrm{d}t} = F_{2}(\mathbf{u}_{i-k},\dots,\mathbf{u}_{i+l}) $$ where $F$ represents a choice of spatial discretization of the two conservation equations, which will depend on neighboring points. The solutions are traveling waves in all four fields. I plan to use a 4th-order Runge-Kutta for the time integration but with varying spatial schemes.
I have read that the stability for nonlinear FD schemes is governed by the local `amplification' matrix. I have no idea how to find the amplification matrix for this system, as the RK4 time step is rather complicated.
Also, can the behavior of the solution be expressed generally in terms of $F$? I expect that the amount of numerical dispersion and diffusion will be related to the accuracy of the spatial scheme, but $F$ by itself is just an expression over the spatial grid points and thus says nothing about how Fourier modes will be affected by the time-stepping scheme.