A scheme for automatic numerical integration is presented, It uses the change of variable x = 1 +(2/'n-){[I + z2)]z%Ii~~_arccos z} to transform the integral to be computed, f~f(x)dx. into (16/37r)f' f(xXl-z2)/ii 7~dz, which is approximated by successive n-points Gauss-Chebyshev quadrature formulas of the second kind (1,,). Due to the special nature of their abscissas and weights, a sequence of formulas I~, 13, 1(~I ~~,, ~+ may be generated, such that '2~+ I may be computed with only n + 1 new integrand evaluations, using the previous value of 1,,. An error estimation is proposed for '2~+I, which only needs two previous values of the sequence (1,, and '(,~+1)/2)The algorithm may be implemented by a very short program (a FORTRAN 77 version is included) that spends practically all its running time in integrand evaluations. It is compared with other methods for automatic numerical integration (trapezoidal rule, Simpson's rule, Romberg's method, an adaptive Gauss-Kronrod rule and Clenshaw-Curtis method) over a broad set of 20 functions. We conclude that the present method is very simple and reliable and is the most efficient among the methods tested here. Possible applications in density functional theory are explored. a poor testing (according to ref.