Historically, the Taylor series representations or truncated Taylor series approximations of a function at a point $x_0$ was first done by taking Newton's form of an interpolation polynomial for points of the form $x_0 + n \Delta$, where $\Delta$ is a positive real number and $n$ is a natural number, and then taking the limit as $\Delta$ goes to $0$. Could someone explain in detail how this was done?
Let $f$ be a function on an interval of real numbers. Let $\Delta$ be a real number and let $x_0$ be in the domain of $f$ such that $x_0$, $x_0 + \Delta$, $\dots$, $x_0 + n\Delta$ is in the domain of $f$, where $n$ is a positive integer. Newton's form of the interpolation polynomial of the data $\{ (x_0 + k \Delta, f(x_0 + k\Delta) \}$ is $$f(x_0) + \frac{f(x_0 + \Delta) - g_0(x_0+\Delta)}{\Delta}(x-x_0) + \cdots + \frac{f(x_0 + n\Delta ) - g_{n-1}(x_0 + n\Delta)}{n!\Delta}(x-x_0)\cdots (x-(n-1)\Delta), $$ where $g_k$ is the interpolation polynomial for the first $k+1$ points. Taking the limit as $\Delta$ tends towards $0$ gives $$ f(x_0) + f'(x_0)\cdot x + \cdots + \frac{1}{n!}\cdot (\lim_{\Delta \to 0} \frac{f(x_0 + n\Delta) -g_{n-1}(x_0 + n\Delta)}{\Delta^n})\cdot x^n, $$ as long as $f$ is sufficently smooth and the limit $\lim_{\Delta \to 0} \frac{f(x_0 + n\Delta) -g_{n-1}(x_0 + n\Delta)}{\Delta^n}$ exits. But why does the limit $\lim_{\Delta \to 0} \frac{f(x_0 + n\Delta) -g_{n-1}(x_0 + n\Delta)}{\Delta^n}$ equal $f^{(n)}(x_0)$?
$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
\begin{align} \mrm{f}\pars{x} & = \mrm{f}\pars{0} + \int_{0}^{x}\mrm{f}'\pars{t}\dd t \,\,\,\stackrel{t\ \mapsto\ x - t}{=}\,\,\, \mrm{f}\pars{0} + \int_{0}^{x}\mrm{f}'\pars{x - t}\dd t \\[5mm] & \stackrel{\mrm{IBP}}{=}\,\,\, \mrm{f}\pars{0} + \mrm{f}'\pars{0}x + \int_{0}^{x}\mrm{f}''\pars{x - t}t\,\dd t \\[5mm] & \stackrel{\mrm{IBP}}{=}\,\,\, \mrm{f}\pars{0} + \mrm{f}'\pars{0}x + {1 \over 2}\,\mrm{f}''\pars{0}x^{2} + {1 \over 2}\int_{0}^{x}\mrm{f}'''\pars{x - t}t^{2}\,\dd t \\[1cm] & \stackrel{\mrm{IBP}}{=}\,\,\, \mrm{f}\pars{0} + \mrm{f}'\pars{0}x + {1 \over 2}\,\mrm{f}''\pars{0}x^{2} + {1 \over 6}\,\mrm{f}'''\pars{0}x^{3} \\[2mm] &\ + {1 \over 6}\int_{0}^{x}\mrm{f}^{\pars{\texttt{IV}}}\pars{x - t}t^{3}\,\dd t \\[1cm] & = \cdots = \bbx{\sum_{k = 0}^{n}\mrm{f}^{\mrm{\pars{k}}}\pars{0}\,{x^{k} \over k!} + {1 \over n!}\ \underbrace{\int_{0}^{x} \mrm{f}^{\mrm{\pars{n + 1}}}\pars{x - t}t^{n}\dd t} _{\ds{\int_{0}^{x} \mrm{f}^{\mrm{\pars{n + 1}}}\pars{t}\pars{x - t}^{n}\,\dd t}}} \end{align}