  This work outlines a fast, high-precision time-domain solver for scalar,
electromagnetic and gravitational perturbations on hyperboloidal foliations of
Kerr space-times. Time-domain Teukolsky equation solvers have typically used
explicit methods, which numerically violate Noether symmetries and are
Courant-limited. These restrictions can limit the performance of explicit
schemes when simulating long-time extreme mass ratio inspirals, expected to
appear in LISA band for 2-5 years. We thus explore symmetric (exponential,
Pad\'e or Hermite) integrators, which are unconditionally stable and known to
preserve certain Noether symmetries and phase-space volume. For linear
hyperbolic equations, these implicit integrators can be cast in explicit form,
making them well-suited for long-time evolution of black hole perturbations.
The 1+1 modal Teukolsky equation is discretized in space using polynomial
collocation methods and reduced to a linear system of ordinary differential
equations, coupled via mode-coupling arrays and discretized (matrix)
differential operators. We use a matricization technique to cast the
mode-coupled system in a form amenable to a method-of-lines framework, which
simplifies numerical implementation and enables efficient parallelization on
CPU and GPU architectures. We test our numerical code by studying late-time
tails of Kerr spacetime perturbations in the sub-extremal and extremal cases.
