The Quasiparticle Random-Phase Approximation (QRPA) based on the nuclear energy density functionals (EDF) is one of the most widely used theoretical frameworks for describing collective excitations in deformed atomic nuclei. Since the configurational space of quasiparticle excitations can become very large, especially for heavy deformed nuclei, the standard matrix eigendecomposition solution of the QRPA equation is often prohibitive in terms of computational resources. Recently, the Finite Amplitude Method (FAM) with its quasiparticle adaptation (QFAM) has been proposed as a feasible method for solving the QRPA equation. Within this doctoral research, a highly efficient implementation of the QFAM solver based on the relativistic nuclear energy density functionals has been developed, namely the DIRQFAM solver in a form of a program package. Due to its efficiency, the developed DIRQFAM solver is suitable for describing collective excitations in axially symmetric atomic nuclei with quadrupole deformed ground states. Even heavy systems and systematic large-scale calculations are within reach. In this dissertation, the implementation of QFAM solver is presented and discussed, together with a new proposed method called the Kernel Polynomial Method (KPM). The KPM method is built on top of QFAM solver and uses QFAM iterations as a means to find the expansion coefficients of the QRPA response function which is expanded in terms of kernel-adjusted Chebyshev polynomials. It is shown that for non-relativistic EDF, the KPM method significantly outperforms the conventional QFAM approach in terms of computational complexity. Two applications of the developed QFAM solver are covered: i) exploration of low-energy response function on multipole excitations in light α-conjugate N = Z nuclei, and ii) description of quasiparticle-vibration coupling in deformed systems where QFAM is used to provide the phonon degrees of freedom for Dyson equation.