Origin of Taylor Series 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)$? 
 A: $\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}$

Besides the historical account ( see @Conifold Pedagogical Answer ), the derivation involves a repeated Integracion by Parts (IBP). Namely,

\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}
