The probability density function (PDF) of the responses of non-linear stochastic system excited by white noise is approximated with the exponential function of polynomial in state variables. Special measure is taken to satisfy FPK equation in the weak sense of integration with the assumed PDF. Gaussian closure method is a special case of the proposed method. Examples are given to show the application of the method to the systems with additive random excitations and those with both additive and multiplicative random excitations. The PDFs obtained with the proposed method and conventional Gaussian closure method are compared with obtainable exact ones. Numerical results showed that the PDFs obtained with the proposed method can be very close to the exact ones regardless of the degree of system non-linearity. In some cases, even exact solution can be obtained with the proposed method.