  Numerical simulations of extereme mass ratio inspirals, the mostimportant
sources for the LISA detector, face several computational challenges. We
present a new approach to evolving partial differential equations occurring in
black hole perturbation theory and calculations of the self-force acting on
point particles orbiting supermassive black holes. Such equations are
distributionally sourced, and standard numerical methods, such as
finite-difference or spectral methods, face difficulties associated with
approximating discontinuous functions. However, in the self-force problem we
typically have access to full a-priori information about the local structure of
the discontinuity at the particle. Using this information, we show that
high-order accuracy can be recovered by adding to the Lagrange interpolation
formula a linear combination of certain jump amplitudes. We construct
discontinuous spatial and temporal discretizations by operating on the
corrected Lagrange formula. In a method-of-lines framework, this provides a
simple and efficient method of solving time-dependent partial differential
equations, without loss of accuracy near moving singularities or
discontinuities. This method is well-suited for the problem of time-domain
reconstruction of the metric perturbation via the Teukolsky or
Regge-Wheeler-Zerilli formalisms. Parallel implementations on modern CPU and
GPU architectures are discussed.
