  The traditional approach to perturbations of nonrotating black holes in
General Relativity uses the reformulation of the equations of motion into a
radial second-order Schr\"odinger-like equation, whose asymptotic solutions are
elementary. Imposing specific boundary conditions at spatial infinity and near
the horizon defines, in particular, the quasi-normal modes of black holes. For
more complicated equations of motion, as encountered for instance in modified
gravity models with different background solutions and/or additional degrees of
freedom, we present a new approach that analyses directly the first-order
differential system in its original form and extracts the asymptotic behaviour
of perturbations, without resorting to a second-order reformulation. As a
pedagogical illustration, we apply this treatment to the perturbations of
Schwarzschild black holes and then show that the standard quasi-normal modes
can be obtained numerically by solving this first-order system with a spectral
method. This new approach paves the way for a generic treatment of the
asymptotic behaviour of black hole perturbations and the identification of
quasi-normal modes in theories of modified gravity.
