  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, such a convenient Schr\"odinger-like reformulation might be
unavailable, even in a generalised matricial form. In order to tackle such
cases, we present a new approach that analyses directly the first-order
differential system in its original form and extracts the asymptotic behaviour
of perturbations. 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.
