Abstract: Probability density evolution method (PDEM) provides a new approach for the dynamic response analysis of nonlinear stochastic structures. One of the key steps in obtaining the probability density function (PDF) of the dynamic response of nonlinear stochastic structures through PDEM is to solve the generalized probability density evolution equation (GDEE). The finite difference method is usually used to solve GDEE. Unfortunately, since GDEE is a first order hyperbolic partial differential equation with variable coefficients and discontinuous initial conditions, solving GDEE by finite difference method may encounter grid sensitivity problems, numerical dispersion and numerical dissipation. In this work, from the perspective of global approximation, a global interpolation scheme was constructed for GDEE based on the Chebyshev pseudo-spectral method, which overcame the problems of numerical dispersion, numerical dissipation, and grid sensitivity. Considering that the coefficients of GDEE remain constant at each time step, a sequential matrix exponential solution for GDEE in the time domain was derived. Since the sequent matrix exponential solution of GDEE is analytically formulated, this method can overcome numerical stability issues. The numerical results from two examples demonstrate that GDEE can be effectively solved using a combination of the sequential matrix exponential method for time integration and the Chebyshev pseudo-spectral method for spatial discretization. The numerical solution exhibited excellent agreement with both the exact solution and the results from Monte Carlo simulations, with negligible numerical dissipation and numerical dispersion. In addition, due to the high convergence of pseudo-spectral method and the fact that the sequential matrix exponential method is not restricted by CFL (Courant-Friedrichs-Lewy) condition, the proposed method has excellent numerical stability and computational efficiency.