Solution 1:

First of all you're right, it is no coincidence. There is a connection between the interpolation formula and the Taylor expansion.

A fine answer to this question can be found in R. Courants Introduction to Calculus and Analysis 1, ch. 5 Taylor's expansion, Appendix.II.3 The Estimate of the Remainder.

In the following I simply take the relevant passages from there together with his nice add-on about a precise meaning to an expression commonly used in geometry.

First step: Notation and preliminary settings which come from Appendix II.2. Construction of the Solution. Newton's Interpolation formula.

There he constructs for a given function $f(x)$ an interpolation polynomial $\phi(x)$ of $n$th degree, such, that $\phi(x_0)=f_0,\dots,\phi(x_n)=f_n$. The construction is done in a stepwise manner starting with an interpolation polynomial $\phi_0(x)=A_0$ with degree zero, then an interpolation polynomial $\phi_1(x)=A_1(x-x_0)$ with degree one which is added to $\phi_0(x)$, till he reaches the polynomial $$\phi(x)=\phi_n(x)=A_0+A_1(x-x_0)+A_2(x-x_0)(x-x_1)+\dots+A_n(x-x_0)\cdots(x-x_{n-1}).$$ The coeffficients $A_i$ can by found by solving the $n+1$ equations \begin{align*} f_0&=A_0\\ f_1&=A_0+A_1(x_1-x_0)\tag{1}\\ \dots&\dots\dots\dots\dots\\ f_n&=A_0+A_1(x_n-x_0)+\dots+A_n(x_n-x_0)(x_n-x_1)\cdots(x_n-x_{n-1})\\ \end{align*}

From R.Courants Introduction to Calculus and Analysis 1, ch. 5 Taylor's expansion, Appendix.II.3: The estimate of the Remainder

In this section $R(x)=f(x)-\phi(x)$, the error of interpolation is introduced. First we only know, that $$R(x_0)=R(x_1)=\cdots=R(x_n)=0$$ In order to be able to say more, $f(x)$ is assumed to have continuous derivatives of at least the $(n+1)$th order.

Now he notes, that for every choice of the constant $c$, the function $$K(x)=R(x)-c(x-x_0)(x-x_1)\cdots(x-x_n)$$ vanishes at the $n+1$ points $x_0,\dots,x_n$. We can then determince $c$ so that $K(y)=0$, that is $$c=\frac{R(y)}{(y-x_0)(y-x_1)\cdots(y-x_n)}$$ So there are $n+2$ points at which $K(x)$ vanishes. Now the generalized Rolle's theorem is applied to $K(x)$; by this we know there is a value $x=\xi$ between the largest and the smallest of the values $x_0,x_1,\dots,x_n,y$, such that $K^{(n+1)}(\xi)=0$. Since $R(x)=f(x)-\phi(x)$, and $\phi$, as a polynomial of $n$th order, has an identically vanishing $(n+1)$th derivative, we have $$f^{(n+1)}(\xi)-c(n+1)!=0$$ noting that $(n+1)!$ is the $(n+1)$th derivative of $(x-x_0)\cdots(x-x_n)$. Thus we have obtained for $c$, a second expression $c=\frac{f^{(n+1)}(\xi)}{(n+1)!}$, containing $\xi$ and depending in some manner on $y$. We now use the equation $K(y)=0$, in which $y$ is completely arbitrary and therefore can be replaced by $x$, and obtain the representation $$R(x)=\frac{(x-x_0)(x-x_1)\cdots(x-x_n)}{(n+1)!}f^{(n+1)}(\xi)$$ where $\xi$ is some value lying between the smallest and the largest of the points $x,x_0,x_1,\dots,x_n$. Thus the general problem of interpolation for a given function $f(x)$ is completely solved. We have for $f(x)$ the representation \begin{align*} f(x)&=A_0+A_1(x-x_0)+A_2(x-x_0)(x-x_1)+\cdots\\ &\qquad\qquad+A_n(x-x_0)(x-x_1)\cdots(x-x_{n-1})+R_n\tag{2} \end{align*} where the coefficients $A_0,A_1,\dots,A_n$ can be found successively from the values of $f$ at the points $x_0,x_1,\dots,x_n$ by the recursion formulas $(1)$ and where the remainder $R_n$ is of the form \begin{align*} R_n=\frac{(x-x_0)(x-x_1)\cdots(x-x_n)}{(n+1)!}f^{(n+1)}(\xi)\tag{3} \end{align*} with a suitable number $\xi$ between the largets and smallest of the values $x,x_0,x_1,\dots,x_n$. If we take the corresponding formula $(2)$ for $f(x)$ with $n$ replaced by $n-1$ and subtract, we obtain $$A_n(x-x_0)(x-x_1)\cdots(x-x_{n-1})+R_n-R_{n-1}=0$$ For $x=x_n$ we have $R_n=0$ and hence for the coefficient $A_n$ (using $(3)$ with $n$ replaced by $n-1$) the representation $$A_n=\frac{f^{(n)}(\xi)}{n!}$$ where $\xi$ lies between the smallest and largest of the values $x_0,x_1,\dots,x_n$. Similar representations exist for $A_{n-1},A_{n-2},\dots,A_0$.

Thus we recognize that if the points $x_0,x_1,\dots,x_n$ are tending together to one and the same point, perhaps the origin, then our interpolation formula (2) goes term for term into the Taylor formula $$f(x)=f(0)+\frac{x}{1!}f'(0)+\frac{x^2}{2!}f''(0)\cdots+\frac{x^n}{n!}f^{(n)}(0)+R_n$$ with the Lagrange form $$R_n=\frac{x^{n+1}}{(n+1)!}f^{(n+1)}(\Theta x)\qquad\qquad0\leq\Theta\leq1$$ of the remainder.

$$ $$

The Taylor formula can thus be considered a limiting case of the Newton interpolation formula.

(And here a nice add-on)

This formula enables us to give precise meaning to an expression commonly used in geometry. The osculating parabola which meets a given curve at a point of $n$th order, is said to have $(n+1)$ consecutive points in common with the given curve at the point. Actually, we obtain this osculating parabola if we find a parabola having $n+1$ points in common with the curve, and then draw these points together.

Analytically, this just corresponds to the transition from the interpolating to the Taylor polynomial. In the same fashion we can characterize the osculation of arbitrary curves. For example, the circle of curvature is that circle which has three consecutive points in common with the given curve.