In this now famous paper by Zabusky and Kruskal, the authors use a finite difference scheme to numerically integrate the Korteweg-de Vries equation:
$$u_t+uu_x+\delta^2u_{xxx}=0$$ The formula they use looks something like (footnote 6): $$u_{n}^{m+1}=u_{n}^{m-1}-\frac{h_t}{3h_x}\left(u_{n+1}^m+u_n^m+u_{n-1}^m\right)\left(u_{n+1}^{m}-u_{n-1}^m\right)-\frac{\delta^2h_t}{h_x^3}\left(u_{n+2}^m-u_{n-2}^m-2\left(u_{n+1}^m-u_{n-1}^m\right)\right)$$ where $h_t,h_x$ are the spacing in time and space, respectively. How did the authors arrive at this formula? I recognize the last term on the RHS as a typical discretization of the third derivative, and it seems like the time discretization is a centered-Euler method, but this doesn't work out to be the same for me.