A new method for numerical integration of a well-behaved function over a finite interval is presented. The method involves expanding the integrand in a series of Chebyshev polynomials and integrating term by term. Examples are given, and the method is compared with Simpson's rule and Gauss' method. The method offers the advantage of easy accuracy checking before integration. It shares some advantages of Simpson's and Gauss' methods, such as the ability to double the number of ordinates without wasting previous work and achieving economy in the number of ordinates. Additionally, it provides values of the indefinite integral throughout the range (a, b), unlike Simpson's rule which only provides values at tabular points, and Gauss' method which is unsuited to indefinite integration. The paper primarily focuses on integrating non-singular functions over finite ranges, but notes that infinite ranges can often be transformed or approximated. Some integrands with weak singularities can also be treated with this method. However, integrals over infinite ranges, such as ∫₀^∞ e^{-x}f(x)dx and ∫_{-∞}^{+∞} e^{-x²}f(x)dx, are best evaluated by special methods like Laguerre-Gauss and Hermite-Gauss formulae. The second integral can also be efficiently evaluated using a simple summation formula. For oscillating functions over infinite ranges, a method based on the Euler transformation of series is described. Simpson's rule is a three-point Newton-Cotes formula, with coefficients that are powers of two. Higher-order formulae have more complex coefficients, which can lead to larger rounding errors. To ensure accuracy, the number of intervals (m) must be sufficiently large. However, for rapidly computable functions, this is often considered worthwhile. For functions with rapidly changing behavior, the interval can be adjusted dynamically based on the agreement of results from different intervals.A new method for numerical integration of a well-behaved function over a finite interval is presented. The method involves expanding the integrand in a series of Chebyshev polynomials and integrating term by term. Examples are given, and the method is compared with Simpson's rule and Gauss' method. The method offers the advantage of easy accuracy checking before integration. It shares some advantages of Simpson's and Gauss' methods, such as the ability to double the number of ordinates without wasting previous work and achieving economy in the number of ordinates. Additionally, it provides values of the indefinite integral throughout the range (a, b), unlike Simpson's rule which only provides values at tabular points, and Gauss' method which is unsuited to indefinite integration. The paper primarily focuses on integrating non-singular functions over finite ranges, but notes that infinite ranges can often be transformed or approximated. Some integrands with weak singularities can also be treated with this method. However, integrals over infinite ranges, such as ∫₀^∞ e^{-x}f(x)dx and ∫_{-∞}^{+∞} e^{-x²}f(x)dx, are best evaluated by special methods like Laguerre-Gauss and Hermite-Gauss formulae. The second integral can also be efficiently evaluated using a simple summation formula. For oscillating functions over infinite ranges, a method based on the Euler transformation of series is described. Simpson's rule is a three-point Newton-Cotes formula, with coefficients that are powers of two. Higher-order formulae have more complex coefficients, which can lead to larger rounding errors. To ensure accuracy, the number of intervals (m) must be sufficiently large. However, for rapidly computable functions, this is often considered worthwhile. For functions with rapidly changing behavior, the interval can be adjusted dynamically based on the agreement of results from different intervals.