The probability density function of the responses of nonlinear random vibration of a multi-degree-of-freedom system is formulated in the defined domain as an exponential function of state variables. The probability density function is assumed to be governed by Fokker-Planck-Kolmogorov (FPK) equation. Special measure is taken to satisfy the FPK equation in the average sense of integration with the assumed function and quadratic algebraic equations are obtained for determining the unknown probability density function. Two-degrees-of-freedom ssytems are analyzed with the proposed method to validate the method for nonlinear multi-degree-of-freedom systems. The probability density functions obtained with proposed method are compared with the obtainable exact and simulated ones. Numerical results showed that the probability density function solutions obtained with the proposed method are much close to the exact and simulated solutions even for highly nonlinear systems with both external and parametric excitations.