  The 3-level leapfrog time integration algorithm is an attractive choice for
numerical relativity simulations since it is time-symmetric and avoids
non-physical damping. In Newtonian problems without velocity dependent forces,
this method enjoys the advantage of long term stability. However, for more
general differential equations, whether ordinary or partial, delayed onset
numerical instabilities can arise and destroy the solution. A known cure for
such instabilities appears to have been overlooked in many application areas.
We give an improved cure ("deloused leapfrog") that both reduces memory demands
(important for 3+1 dimensional wave equations) and allows for the use of
adaptive timesteps without a loss in accuracy. We show both that the
instability arises and that the cure we propose works in highly relativistic
problems such as tightly bound geodesics, spatially homogeneous spacetimes, and
strong gravitational waves. In the gravitational wave test case (polarized
waves in a Gowdy spacetime) the deloused leapfrog method was five to eight
times less CPU costly at various accuracies than the implicit Crank-Nicholson
method, which is not subject to this instability.
