A bivariate Chebyshev polynomials approach is proposed to estimate the dynamic response bounds of nonlinear systems with interval uncertainties. The existing collocation method directly searches the maximum and minimum values of the surrogate model in the entire interval space by the scanning method (SM). The presence of too many uncertain parameters will lead to expansive computational cost. To overcome this shortcoming, the dynamic response is decomposed by a bivariate function decomposition (BFD), established based on high-order Taylor expansion, into the sum of multiple univariate and bivariate response functions. The above univariate and bivariate functions are fitted using Chebyshev polynomials, and polynomial coefficients are obtained through one-dimensional (1D) and two-dimensional (2D) interpolation points. Thus, the solution of the nonlinear dynamic systems with uncertain parameters can be transformed into that of univariate and bivariate Chebyshev interval functions. The extremum values of the low-dimensional Chebyshev interval functions can be found by SM, and then the bounds of dynamic response are acquired by interval arithmetic. Since SM searches for extreme values only in 1D and 2D uncertain domains, the amount of calculation is reduced compared to searching the whole uncertain space. The efficiency, practicability and effectiveness of the proposed interval uncertainty analysis method are proved by three dynamic examples.