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.