We present a numerical approach to calculate non-equilibrium eigenstates of a periodically time-modulated quantum system. The approach is based on the use of a chain of single-step time-independent propagating operators. Each operator is time-specific and constructed by combining the Magnus expansion of the time-dependent system Hamiltonian with the Chebyshev expansion of an operator exponent. A construction of a unitary matrix of the Floquet operator, which evolves a system state over the full modulation period, is performed by propagating the identity matrix over the period. The independence of the evolutions of basis vectors makes the propagation stage suitable for implementation on a parallel cluster. Once the propagation stage is completed, a routine diagonalization of the Floquet matrix is performed. Finally, an additional propagation round, now with the eigenvectors as the initial states, allows to resolve the time-dependence of the Floquet states and calculate their characteristics. We demonstrate the accuracy and scalability of the algorithm by applying it to calculate the Floquet states of two quantum models, namely (i) a synthesized random-matrix Hamiltonian and (ii) a many-body Bose-Hubbard dimer, both of the size up to 10 4 states.the Hamiltonian H(t) but instead of the unitary Floquet operator