This paper presents the development of an efficient and numerically stable algorithm for accurately computing the eigensolutions of the quadratic real symmetric eigenproblem arising in the finite dynamic element (FDE) formulation. Closely related to the subspace forward iteration method, the proposed scheme is well suited to extracting the lowest natural frequencies and associated mode shapes of large practical eigenproblems, takes full advantage of the banded configuration of the stiffness, inertia and dynamic correction matrices involved in the eigenproblem and ensures a monotonic convergence from above to the required eigenpairs. A shifting technique for convergence acceleration and an eigenpair verification scheme are also presented. Numerical examples are shown demonstrating the excellent performance of the solution procedure.