Add correct initialization to second order differencing scheme: without this, the norm is initially oscillating, unless one goes to ridiculously short time steps. So far (Jan 2023), this has escaped our attention over the years because none of the demo examples uses this integrator.
Add 4th order schemes, in analogy to the WaveTrain codes where this scheme has been shown to be very good.
Don't forget to explain this method here