The paper proposes a new methodology for time-dependent reliability analysis of vibratory systems using a combination of a first-order, four-moment (FOFM) method and a non-Gaussian Karhunen–Loeve (NG-KL) expansion. The approach can also be used for random vibrations studies. The vibratory system is nonlinear and is excited by stationary non-Gaussian input random processes which are characterized by their first four marginal moments and autocorrelation function. The NG-KL expansion expresses each input non-Gaussian process as a linear combination of uncorrelated, non-Gaussian random variables and computes their first four moments. The FOFM method then uses the moments of the NG-KL variables to calculate the moments and autocorrelation function of the output processes based on a first-order Taylor expansion (linearization) of the system equations of motion. Using the output moments and autocorrelation function, another NG-KL expansion expresses the output processes in terms of uncorrelated non-Gaussian variables in the time domain, allowing the generation of output trajectories. The latter are used to estimate the time-dependent probability of failure using Monte Carlo simulation (MCS). The computational cost of the proposed approach is proportional to the number of NG-KL random variables and is significantly lower than that of other recently developed methodologies which are based on sampling. The accuracy and efficiency of the proposed methodology is demonstrated using a two-degree-of-freedom nonlinear vibratory system with random coefficients excited by a stationary non-Gaussian random process.