This paper presents a new approach for calculation of uncertainty in dynamic simulation results. The statistical moments (mean, variance, skewness etc.) of the simulation results are calculated using Gaussian-quadrature with "customized" weight function. Based on these moments, an approximating probability density function (pdf) is created by expansion into orthogonal polynomial series. The percentiles of the distribution can then be calculated. The method is computationally less demanding than Monte-Carlo simulation when the number of uncertain parameters are limited. A number of examples are used to illustrate the applicability of the proposed framework.