Estimate Taylor series coefficients from samples of a function
Instead of exact polynomial fitting, you could use a least-squares fit, which will find the polynomial of specified order that minimizes the total squared error between the fit and the measured $(x_i,y_i)$ pairs. This can help to mitigate the effects of noise on the fit.

Given measurements $y_i$ of a function $y = f(x)$ at domain values $x_i$ ($i = 0, 1, \ldots , N$), choose a polynomial order $M \le N$ (if $M = N$, then you're down to exact polynomial fitting, as the $N$ points uniquely determine an $M$th order polynomial). Then, set up a system of equations that are linear in the desired polynomial coefficients $p_k$:

$$
y_i = p_Mx_i^M + p_{M-1}x_i^{M-1} + \ldots + p_1x_i + p_0, \;\;\;\;\;i = 0, 1, \ldots , N
$$

The least-squares problem can be solved by arranging the measurements into matrix-vector form:

$$
\mathbf{A} =  \left [ {\begin{array}{cc}
 x_0^M &amp; x_0^{M-1} &amp; \cdots &amp; x_0 &amp; 1  \\
 x_1^M &amp; x_1^{M-1} &amp; \cdots &amp; x_1 &amp; 1  \\
 \vdots &amp; \vdots &amp; \vdots &amp; \vdots &amp; \vdots \\
 x_N^M &amp; x_N^{M-1} &amp; \cdots &amp; x_N &amp; 1 \end{array} } \right ], 
\;\;\mathbf{y}= \left [ {\begin{array}{cc}
 y_0  \\
 y_1 \\
 \vdots  \\
 y_N  \end{array} } \right ]
$$

The least-squares solution generates the vector of polynomial coefficients $\left[p_M, p_{M-1}, \ldots , p_0\right]$ that minimizes the total squared error in the above linear system. The solution can be computed as:

$$
\mathbf{\tilde p} = (\mathbf{A^TA})^{-1} \mathbf{A^Ty}
$$

It is worth noting that the matrix $(\mathbf{A^TA})^{-1} \mathbf{A^T}$ is also known as the pseudoinverse of the matrix $\mathbf{A}$. You can then use the least-squares polynomial coefficient vector $\mathbf{\tilde p}$ to evaluate the polynomial at any other $x$ values that you wish.