The main culprit is the low variability of the dependent variable: `tableofvalues[[All,5]]`. If you standardize that variable (subtracting the mean and then dividing by the standard deviation), then all your troubles go away. In fact standardizing variables is almost always recommended even for the independent variables whenever there might be round-off error involved especially for the fitting of polynomials. However, the suggestion by @MarcoB (setting a higher precision) achieves essentially the same result. stdvalues = tableofvalues; stdvalues[[All, 5]] = Standardize[tableofvalues[[All, 5]]]; stdev5 = StandardDeviation[tableofvalues[[All, 5]]]; polynomtabelle = {1, x[1], x[1]^2, x[1]^3, x[1]^4, x[1]^5, x[2], x[2]^2, x[2]^3, x[2]^4, x[2]^5, x[3], x[3]^2, x[3]^3, x[3]^4, x[3]^5, x[4], x[4]^2, x[4]^3, x[4]^4, x[4]^5} polynom = LinearModelFit[stdvalues, polynomtabelle, Array[x, {4}]]; NumberForm[ stdev5 Extract[ CoefficientList[polynom["BestFit"], Array[x, {4}]], {5, 0, 0, 0} + 1], Floor@$MachinePrecision] NumberForm[ stdev5 Extract[ CoefficientList[polynom["BestFit"], Array[x, {4}]], {0, 0, 0, 5} + 1], Floor@$MachinePrecision] with output (* 0.0619247172791471 *) (* 0.0619247172736994 *) Note that one must multiply by the standard deviation of `tableofvalues[[All,5]]` to get the appropriate coefficient.