  Gravitational wave signals from extreme mass ratio inspirals are a key target
for space-based gravitational wave detectors. These systems are typically
modeled as a distributionally-forced Teukolsky equation, where the smaller
black hole is treated as a Dirac delta distribution. Time-domain solvers often
use regularization approaches that approximate the Dirac distribution that
often introduce small length scales and are a source of systematic error,
especially near the smaller black hole. We describe a multi-domain
discontinuous Galerkin method for solving the distributionally-forced Teukolsky
equation that describes scalar fields evolving on a Kerr spacetime. To handle
the Dirac delta, we expand the solution in spherical harmonics and recast the
sourced Teukolsky equation as a first-order, one-dimensional symmetric
hyperbolic system. This allows us to derive the method's numerical flux to
correctly account for the Dirac delta. As a result, our method achieves global
spectral accuracy even at the source's location. To connect the near field to
future null infinity, we use the hyperboloidal layer method, allowing us to
supply outer boundary conditions and providing direct access to the far-field
waveform. We document several numerical experiments where we test our method,
including convergence tests against exact solutions, energy luminosities for
circular orbits, the scheme's superconvergence properties at future null
infinity, and the late-time tail behavior of the scalar field. We also compare
two systems that arise from different choices of the first-order reduction
variables, finding that certain choices are numerically problematic in
practice. The methods developed here may be beneficial when computing
gravitational self-force effects, where the regularization procedure has been
developed for the spherical harmonic modes and high accuracy is needed at the
Dirac delta's location.
