The quasinormal modes of relativistic compact objects encode important information about the gravitational response associated with astrophysical phenomena. Detecting such oscillations would provide us with a unique understanding of the properties of compact stars and may give definitive evidence for the existence of black holes. However, computing quasinormal modes in realistic astrophysical environments is challenging due to the complexity of the spacetime background and of the dynamics of the perturbations. We discuss two complementary methods for computing the quasinormal modes of spherically symmetric astrophysical systems, namely, the direct integration method and the continued-fraction method. We extend these techniques to dealing with generic coupled systems of linear equations, with the only assumption being that the interaction between different fields is effectively localized within a finite region. In particular, we adapt the continued-fraction method to include cases where a series solution can be obtained only outside an effective region. As an application, we compute the polar quasinormal modes of boson stars by using the continued-fraction method for the first time. The methods discussed here can be applied to other situations in which the perturbations effectively couple only within a finite region of space.