This paper presents a numerical integration method for estimating the area under the curve (AUC) over the infinite time interval. This method is based on the Gauss-Laguerre quadrature and produces AUC estimates over the infinite time interval without extrapolation in a usual sense. By contrast, in traditional schemes, piecewise interpolation is used to obtain the area up to the final sampling point, and the remaining portion is extrapolated using nonlinear regression. In this case, there is no theoretical consistency between the quadrature and extrapolation. The inconsistency may cause certain problems. For example, the optimal sampling criterion for the former is not necessarily optimal for the latter. Such inconsistency does not arise in the method of this work. The sampling points are placed near the zeros of Laguerre polynomials so as to directly estimate the AUC over the infinite time interval. The sampling design requires no particular prior information. This is also advantageous over the previous strategy, which worked by minimizing the variance of estimated AUC under the assumptions of particular pharmacokinetic and variance functions. The original Gaussian quadrature is belived to be inappropriate for numerical integration of data because of several restrictions. In this paper, it is shown that, using a simple strategy for managing errors due to these restrictions, the method produces an estimate of AUC with practically sufficient precision. The efficacy of this method is finally shown by numerical simulations in which the bias and variance of its estimate were compared with those of the previous methods such as the trapezoidal, log-trapezoidal, Lagrange, and parabolas-through-the-origin methods.