  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
non-spherical geometries. In this work, 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 $^{38}$Si and $^{250}$Cf are presented and discussed.
