This paper reviews various approximation methods for evaluating one form of the exponential integral known as the Theis well function W( u) and presents an algorithm that can be easily implemented to one's desired accuracy for 0 < u < ∞. The algorithm is based on a combination of a fast-converging series representation for small u and an easy-implementing Gauss-Laguerre quadrature formula when u becomes large. The partition point for u varies with the desired program accuracy and efficiency. The error analyses suggested that a partition point in the interval of 1 ≤ u ≤ 4 can yield sufficiently accurate results for practical applications of groundwater problems with minimum computational effort. A comparison of the results revealed that the popular Stehfest inverse Laplace transform technique generally yielded a relative error several orders of magnitude greater than that obtained by the proposed algorithm. For large u, the approximation accuracy increased rapidly using the quadrature formula while the results of the Laplace inversion technique showed numerical oscillation between some positive and negative values.