We systematically investigate the quasinormal modes of thick branes in $\textit{f}$($\textit{R}$) gravity by numerically solving the Schrödinger-like perturbation equation of gravitational perturbations. To ensure reliability of the results, we employ three complementary methods: the asymptotic iteration method, direct integration of the wave equation, and time-domain numerical evolution. We analyze how model parameters influence the shape of the effective potential of gravitational perturbations and find that the structure of the potential barrier plays a significant role in shaping the quasinormal frequency spectrum. The results obtained from the three methods exhibit strong consistency, thereby ensuring the reliability of the calculations. In particular, real parts of the quasinormal frequencies exhibit an approximately arithmetic progression, suggesting that quasi-localized states can be understood as resonances between the barriers.