This paper proposes a non-intrusive computational method for mechanical dynamic systems involving a large-scale of interval uncertain parameters, aiming to reduce the computational costs and improve accuracy in determining bounds of system response. The screening method is firstly used to reduce the scale of active uncertain parameters. The sequential high-order polynomials surrogate models are then used to approximate the dynamic system’s response at each time step. To reduce the sampling cost of constructing surrogate model, the interaction effect among uncertain parameters is gradually added to the surrogate model by sequentially incorporating samples from a candidate set, which is composed of vertices and inner grid points. Finally, the points that may produce the bounds of the system response at each time step are searched using the surrogate models. The optimization algorithm is used to locate extreme points, which contribute to determining the inner points producing system response bounds. Additionally, all vertices are also checked using the surrogate models. A vehicle nonlinear dynamic model with 72 uncertain parameters is presented to demonstrate the accuracy and efficiency of the proposed uncertain computational method.