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 following the traditional translation process model. 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 cumulative distribution function (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 EOMs 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 are demonstrated using a nonlinear Duffing oscillator problem under a quadratic random excitation and a nonlinear half-car suspension example.