Starting from a general many-body fermionic Hamiltonian, we derive the equations of motion (EOM) for nucleonic propagators in a superfluid system. The resulting EOM is of the Dyson type formulated in the basis of Bogoliubov's quasiparticles. As the leading contributions to the dynamical kernel of this EOM in strongly- coupled regimes contain phonon degrees of freedom in various channels, an efficient method of calculating phonon's characteristics is required to successfully model these kernels. The traditional quasiparticle random phase approximation (QRPA) solvers are typically used for this purpose in nuclear structure calculations, however, they become very prohibitive in nonspherical geometries. In this paper, by linking the notion of the quasiparticle- phonon vertex to the variation of the Bogoliubov's Hamiltonian, we show that the recently developed finite-amplitude method (FAM) can be efficiently employed to compute the vertices within the FAM- QRPA. To illustrate the validity of the method, calculations based on the relativistic density- dependent point-coupling Lagrangian are performed for the single-nucleon states in heavy and medium- mass nuclei with axial deformations. The cases of 38Si and 250Cf are presented and discussed.