The probability density function (PDF) solution of the response is formulated for nonlinear systems under discrete Poisson impulse excitation. The PDF solution is governed by the Kolmogorov–Feller (KF) equation, which is approximately solved by the exponential–polynomial closure (EPC) method. A Duffing oscillator is further investigated in the case of either Gaussian or non-Gaussian distributed amplitude of Poisson impulse to show the effectiveness of the EPC method in these cases. The numerical analysis shows that the EPC method with the polynomial order being 6 presents a good result compared with the simulated result, even in the tails of the PDF of the oscillator response.