This paper proposes a new computationally efficient methodology for random vibrations of nonlinear vibratory systems using a time-dependent second-order adjoint variable (AV2) method, and a second-order projected differentiation (PD2) method. The proposed approach is called AV2-PD2. The vibratory system can be excited by stationary Gaussian or non-Gaussian random processes. A Karhunen-Loeve (KL) expansion expresses each input random process in terms of standard normal random variables. A second-order adjoint approach is used to obtain the required first and second-order output derivatives accurately by solving as many sets of equations of motion (EOMs) as the number of KL random variables. These derivatives are used to compute the marginal CDF of the output process with second-order accuracy. Then, a second-order projected differentiation method calculates the autocorrelation function of each output process with second-order accuracy, at an additional cost of solving as many sets of EOM as the number of outputs of interest, independently of the time horizon (simulation time). The total number of solutions of the EOM scales linearly with the number of input KL random variables and the number of output processes. The efficiency and accuracy of the proposed approach is demonstrated using a non-linear Duffing oscillator problem under a quadratic random excitation.