Numerical methods for evaluating higher-order quantum-mechanical phase integrals are developed. The integrand of each contour integral is factorized inside the contour into an analytic and a singular part, the combination then being treated in the manner of product integration. After transformation the analytic part F(z) is approximated on a standard interval by a finite Chebyshev series sum and the coefficients are evaluated by Gauss-Chebyshev or Lobatto-Chebyshev quadrature (equivalent to Lagrangian interpolation at the Chebyshev or the Clenshaw-Curtis points). Themodified moments or weight integrals giving the contributions to the contour integral from individual Chebyshev polynomial components of F(z) are evaluated by analytically derived recursions. The methods are applied to a quantum scattering model in physical chemistry.