  A hybrid method is developed based on the spectral and finite-difference
methods for solving the inhomogeneous Zerilli equation in time-domain. The
developed hybrid method decomposes the domain into the spectral and
finite-difference domains. The singular source term is located in the spectral
domain while the solution in the region without the singular term is
approximated by the higher-order finite-difference method.
  The spectral domain is also split into multi-domains and the
finite-difference domain is placed as the boundary domain. Due to the global
nature of the spectral method, a multi-domain method composed of the spectral
domains only does not yield the proper power-law decay unless the range of the
computational domain is large. The finite-difference domain helps reduce
boundary effects due to the truncation of the computational domain. The
multi-domain approach with the finite-difference boundary domain method reduces
the computational costs significantly and also yields the proper power-law
decay.
  Stable and accurate interface conditions between the finite-difference and
spectral domains and the spectral and spectral domains are derived. For the
singular source term, we use both the Gaussian model with various values of
full width at half maximum and a localized discrete $\delta$-function. The
discrete $\delta$-function was generalized to adopt the Gauss-Lobatto
collocation points of the spectral domain.
  The gravitational waveforms are measured. Numerical results show that the
developed hybrid method accurately yields the quasi-normal modes and the
power-law decay profile. The numerical results also show that the power-law
decay profile is less sensitive to the shape of the regularized
$\delta$-function for the Gaussian model than expected. The Gaussian model also
yields better results than the localized discrete $\delta$-function.
