Based on the theory of Absolute Nodal Coordinate Formulation (ANCF), this paper proposes a new dynamic computation method to solve the flexible multibody system with uncertain material properties (Young's modulus and Poisson's ratio) that may be induced by the material asymmetric distribution.Rather than traditionally considering an uncertain factor as one single variable in the whole system, the material properties vary continuously in the space domain so that they are described by the random field, which is then discretized to countable random variables using the expansion optimization linear estimation (EOLE) method. The uncertain response of the system is approximated by the Polynomial Chaos (PC) expansion, numerically implemented by a collocation method. We propose and prove an important theory that the collocation method provides the same results as the Gaussian quadrature formula if the roots of the corresponding orthogonal polynomials are used as the collocation points. As a result, the proposed method is a non-intrusive technique that does not modify the original solver but only adds a pre-process and a post-process. The uncertain displacement of system is finally illustrated by an ellipse and ellipsoid, which visually show the uncertainty extent and correlation between different coordinates. The numerical examples show that the proposed method has almost an equivalent accuracy of Monte Carlo simulation but much higher efficiency.