The advent of computing has enabled scientists, engineers, and designers to use numerical methods to calculate an approximate solution of the equations that describe vibratory phenomena. The Finite Element Method is among the most applied tools to perform the harmonic analysis of structures. However, the resulting models may become large in terms of the total number of degrees of freedom, especially when they are highly complex and aiming for high precision solutions. The evaluation of typical models for vibration analysis with sufficient accuracy for engineering purposes can take weeks in modern computing systems, which often restricts their practical applications. It motivates the development of model order reduction techniques as a means of decreasing the solution processing time. The computational cost is reduced by dramatically diminishing the number of degrees of freedom of the numerical model without significantly compromising the solution’s accuracy. The main idea is to capture the most relevant information from the dynamic system with fewer degrees of freedom as possible. This work presents the model order reduction method based on projection in second-order Krylov subspaces and proposes an automatic reduction procedure, which reduces the user intervention and guarantees solution quality by error estimators. The method is implemented in an open-source code using Python. A beam model and a plate model, both under point load, are used for validation. The speed-up factor, regarding the direct solution method, reached hundreds of times for both models. The proposed method is suitable for systems analyzed in a wide frequency range, whose processing time for the full order model solution is high.