  Partial differential equations exhibiting an anisotropic scaling between
space and time -- such as those of Horava-Lifshitz gravity -- have a dispersive
nature. They contain higher-order spatial derivatives, but remain second order
in time. This is inconvenient for performing long-time numerical evolutions, as
standard explicit schemes fail to maintain convergence unless the time step is
chosen to be very small. In this work, we develop an implicit evolution scheme
that does not suffer from this drawback, and which is stable and second-order
accurate. As a proof of concept, we study the numerical evolution of a Lifshitz
scalar field on top of a spherically symmetric black hole space-time. We
explore the evolution of a static pulse and an (approximately) ingoing
wave-packet for different strengths of the Lorentz-breaking terms, accounting
also for the effect of the angular momentum eigenvalue and the resulting
effective centrifugal barrier. Our results indicate that the dispersive terms
produce a cascade of modes that accumulate in the region in between the Killing
and universal horizons, indicating a possible instability of the latter.
