Implementation [of my PDEs]

All the numerical solutions shown were found using the NDSolve function built into Mathematica. In general, finite difference methods, the method of lines and pseudospectral methods can be used. For equations of the form

∂_ttu[t, x]  ∂_xxu[t, x] + f[u[t, x]]

one can set up a simple finite difference method by taking f in the form of pure function and creating from it a kernel with space step dx and time step dt:

PDEKernel[f_, {dx_, dt_}] := Compile[{a,b,c,d}, Evaluate[(2 b - d) + ((a + c - 2 b)/dx² + f[b]) dt²]]

Iteration for n steps is then performed by

PDEEvolveList[ker_, {u0_, u1_}, n_] := Map[First, NestList[PDEStep[ker, #]&, {u0, u1}, n]]

PDEStep[ker_, {u1_, u2_}] := {u2, Apply[ker, Transpose[{RotateLeft[u2], u2, RotateRight[u2], u1}], {1}]}

With this approach an approximation to the top example on page 165 can be obtained from

PDEEvolveList[PDEKernel[(1 - #²)(1 + #)&, {.1, .05}], Transpose[Table[{1, 1} N[Exp[-x²]], {x, -20, 20, .1}]], 400]

For both this example and the middle one the results converge rapidly as dx decreases. But for the bottom example, the pictures below show that convergence is not so rapid, and indeed, as is typical in working with PDEs, despite having used large amounts of computer time I do not know whether the details of the picture in the main text are really correct. The energy function (see above) is at least roughly conserved, but it seems quite likely that the "shocks" visible are merely a consequence of the discretization procedure used.