We present an efficient method for calculating strength functions using the finite amplitude method (FAM) for deformed superfluid heavy nuclei within the framework of the nuclear density functional theory. We demonstrate that FAM reproduces strength functions obtained with the fully self-consistent quasi-particle random-phase approximation (QRPA) at a fraction of computational cost. As a demonstration, we compute the isoscalar and isovector monopole strength for strongly deformed configurations in 100 Zr and 240 Pu by considering huge quasi-particle QRPA spaces. Our approach to FAM, based on Broyden's iterative procedure, opens the possibility for large-scale calculations of strength distributions in well-deformed superfluid nuclei across the nuclear landscape.