When considering the propagation of waves on shallow water surfaces, a simple model is the KdV (Korteweg de Vries) equation :
$$\zeta_t + \zeta_x + \frac{3 \varepsilon}{2} \zeta \zeta_x + \frac{\varepsilon}{6} \zeta_{xxx} = 0$$
where $\zeta_x$ denotes the derivative of $\zeta$ with respect to $x$. $\zeta(x,t)$ is the scaled elevation of the surface of the wave and is a function of the time $t$ and the space $x$. The true elevation is equal to $h \times (1 + \varepsilon \zeta)$, the parameter $\varepsilon$ will be chosen less than $1$ (e.g. $0.1$). We will reduce the computations on the interval $[-L, L]$ and impose periodic boundary conditions.
Let $U$ be a vector of size $N$, discretizing $\zeta$ on a regular subdivision $x_k=\frac{2k-N+1}{N-1} L$ of the interval $[-L, L]$ ($U_k = \zeta(x_k, t)$). We are going to compute a sequence of vector $U^n$ representing the temporal evolution of the wave. The initial condition will be set to a gaussian function : $\zeta(x, 0) = \exp(-\alpha x^2)$ with $\alpha > 0$. As for the heat equation, a finite difference scheme is used to solve this equation. In the following, we describe an algorithm to compute $U^{n+1}$ as a function of $U^n$.
Periodic boundary conditions induce the following conditions :
$$\begin{cases} U_{-1} = U_{N-1} \\\ U_{-2} = U_{N-2} \\\ U_{N} = U_{0} \\\ U_{N+1} = U_{1} \end{cases}$$
We denote $G$ the application which transforms the vector $U$ into the vector :
$$[ G(U) ]_i = \frac{U_{i+1} - U_{i-1}}{2 \Delta x} + \frac{\varepsilon}{4 \Delta x} (U_{i+1} - U_{i-1})(U_{i-1} + U_i + U_{i+1}) + \frac{\varepsilon}{12 \Delta x^3} (U_{i+2} - 2 U_{i+1} + 2 U_{i-1} - U_{i-2})$$
This application is non-linear and can not be represented as a matrix. We have the following evolution system : $\frac{dU}{dt} + G(U) = 0$. In order to solve this evolution equation, we consider the following equation :
$$\frac{U^{n+1} - U^n}{\Delta t} + G\Big(\frac{U^n + U^{n+1}}{2}\Big) = 0$$