  Black hole solutions in general relativity are simple. The frequency spectrum
of linear perturbations around these solutions (i.e., the quasinormal modes) is
also simple, and therefore it is a prime target for fundamental tests of black
hole spacetimes and of the underlying theory of gravity. The following
technical calculations must be performed to understand the imprints of any
modified gravity theory on the spectrum: 1. Identify a healthy theory; 2. Find
black hole solutions within the theory; 3. Compute the equations governing
linearized perturbations around the black hole spacetime; 4. Solve these
equations to compute the characteristic quasinormal modes. In this work (the
first of a series) we assume that the background spacetime has spherical
symmetry, that the relevant physics is always close to general relativity, and
that there is no coupling between the perturbation equations. Under these
assumptions, we provide the general numerical solution to step 4. We provide
publicly available data files such that the quasinormal modes of {\em any}
spherically symmetric spacetime can be computed (in principle) to arbitrary
precision once the linearized perturbation equations are known. We show that
the isospectrality between the even- and odd-parity quasinormal mode spectra is
fragile, and we identify the necessary conditions to preserve it. Finally, we
point out that new modes can appear in the spectrum even in setups that are
perturbatively close to general relativity.
