Random vibration analysis is a mathematical tool that offers great advantages in predicting the mechanical response of structural systems subjected to external dynamic loads whose nature is intrinsically stochastic, as in cases of sea waves, wind pressure, and vibrations due to road asperity. Using random vibration analysis is possible, when the input is properly modeled as a stochastic process, to derive pieces of information about the structural response with a high quality (if compared with other tools), especially in terms of reliability prevision. Moreover, the random vibration approach is quite complex in cases of non-linearity cases, as well as for non-stationary inputs, as in cases of seismic events. For non-stationary inputs, the assessment of second-order spectral moments requires resolving the Lyapunov matrix differential equation. In this research, a numerical procedure is proposed, providing an expression of response in the state-space that, to our best knowledge, has not yet been presented in the literature, by using a formal justification in accordance with earthquake input modeled as a modulated white noise with evolutive parameters. The computational efforts are reduced by considering the symmetry feature of the covariance matrix. The adopted approach is applied to analyze a multi-story building, aiming to determine the reliability related to the maximum inter-story displacement surpassing a specified acceptable threshold. The building is presumed to experience seismic input characterized by a non-stationary process in both amplitude and frequency, utilizing a general Kanai–Tajimi earthquake input stationary model. The adopted case study is modeled in the form of a multi-degree-of-freedom plane shear frame system.