Response of nonlinear systems subjected to harmonic, parametric, and random excitations is of importance in the field of structural dynamics. The transitional probability density function (PDF) of the random response of nonlinear systems under white or colored noise excitation (delta correlated) is governed by both the forward Fokker–Planck (FP) and the backward Kolmogorov equations. This paper presents a new approach for efficient numerical implementation of the path integral (PI) method in the solution of the FP equation for some nonlinear systems subjected to white noise, parametric, and combined harmonic and white noise excitations. The modified PI method is based on a non-Gaussian transition PDF and the Gauss–Legendre integration scheme. The effects of white noise intensity, amplitude, and frequency of harmonic excitation and the level of nonlinearity on stochastic jump and bifurcation behaviors of a hardening Duffing oscillator are also investigated.