General relativity can be tested at many scales using various physical systems. A particularly interesting probe is the study of the ringdown phase of a binary black hole merger, during which the newly-formed black hole emits gravitational waves at given frequencies called its quasinormal modes. Such modes depend heavily on the theory of gravity underlying the solution and can thus be used to test GR and put constraints on modified gravity theories. Quasinormal modes are computed perturbatively: one must obtain the equation of propagation of each degree of freedom of gravitational waves, and this can prove quite difficult in theories where these degrees couple together (scalar-tensor theories for example). We present a systematic approach that extracts the asymptotic behaviour of perturbations at spatial infinity and near the horizon directly from the perturbed Einstein’s equations. This asymptotic behavior gives physical insight about the modes, for example their direction of propagation, and allows us to numerically compute their frequencies, as we show for a particular nonrotating black hole solution of Horndeski’s theory of gravity.