<p>We present a method for numerical modeling of the properties of linear active oscillation modes in gyrotrons, i.e., powerful vacuum generators of electromagnetic waves, which operate in the range from ultra-high to terahertz frequencies. The method is based on the mathematical model of the linear stationary gyrotron theory within the approximation of an unfixed field structure. To formulate a linear spectral problem, we use the method of linearization of the boundary conditions of radiation over the frequency. Using the finite difference method, the boundary value problem is reduced to a linear generalized matrix eigenvalue problem, which is solved by the shift-invert Arnoldi method. An algorithm is proposed for sequentially calculating the spectrum of complex frequencies of linear modes in a gyrotron to find a given number of modes without an additional procedure of searching for initial approximations. The developed technique ensures the calculation of the properties of a given number of linear oscillation modes in a gyrotron, specifically, frequencies, increments (decrements) of instabilities and distributions of the complex field amplitudes, and the current grouped in the electron beam along the cavity. Additionally, the instability boundaries on the plane of the corresponding parameters (starting generation conditions) are calculated. The data of test calculations for an actual gyrotron design are presented. The proposed methods and algorithms will significantly facilitate the process of developing gyrotrons intended for various practical applications.</p>