<p>This paper proposes a Cholesky-parameterized path integration method for numerically solving the transient multi-peak probability density function (PDF) of nonlinear dynamical systems. The proposed path integration method discretizes the transient PDF in the spatial and temporal domains. The Gauss-Legendre quadrature scheme and short-time Gaussian approximation technique fulfill the solution task. The Gaussian transition PDF between two successive time steps is obtained by the Cholesky-parameterized moment equations. This formulation ensures that the covariance matrix keeps positive definiteness for the transition Gaussian PDF when strong excitation or strong nonlinearity exists. The solution procedure is applied to tri-stable Duffing oscillators subjected to three types of Gaussian white noise excitations: purely external, combined external and parametric, and correlated excitations. Both stationary and nonstationary PDFs are obtained at representative time instants. Comparisons with exact stationary solutions and Monte Carlo simulations demonstrate excellent accuracy and significantly reduced computational cost. The numerical analysis also reveals distinct PDF evolution patterns: under purely external excitation, a symmetric triple-peak PDF emerges with a dominant central peak. When an independent parametric excitation is introduced, the symmetry is preserved, while the peak disparity is accentuated. When excitations are correlated, symmetry breaks and asymmetry is developed in the peak distribution. The framework provides an efficient and reliable computational tool for capturing transient multi-peak behavior and offers insights into noise-induced phenomena in multi-stable systems.</p>