[en] A new iterative procedure in fitting data is proposed which can properly correct small nonlinear effects. Applying this procedure to fit the Polyakov loop correlation functions previously obtained in lattice QCD simulations, we are able to eliminate the finite-size boundary effects and obtain the quark potential functions themselves, instead of just parameters for the potentials. The quality of the fits is improved significantly. The potential functions obtained at β=6.0--6.2 scale in fine detail, providing important further support of the scaling. The improved estimates of the string tension satisfy the asymptotic scaling within 97%. In comparison with experiments, the obtained value √σ /Λ_L=79±3 indicates Λ_MS=449±17 MeV, where MS denotes the minimal subtraction scheme. The coupling constant α=0.43±0.03 agrees very well with experiments, and is also in excellent agreement with the theoretical value π/8=0.39, assuming that the long-distance behavior of QCD can be effectively described by the Neveu-Schwarz string theory