Csdrhrt

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 + nDelta$ 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 + kDelta) }$ is
$$f(x_0) + frac{f(x_0 + Delta) - g_0(x_0+Delta)}{Delta}(x-x_0) + cdots + frac{f(x_0 + nDelta ) - g_{n-1}(x_0 + nDelta)}{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 + nDelta) -g_{n-1}(x_0 + nDelta)}{Delta^n})cdot x^n,
$$
as long as $f$ is sufficently smooth and the limit $lim_{Delta to 0} frac{f(x_0 + nDelta) -g_{n-1}(x_0 + nDelta)}{Delta^n}$ exits. But why does the limit $lim_{Delta to 0} frac{f(x_0 + nDelta) -g_{n-1}(x_0 + nDelta)}{Delta^n}$ equal $f^{(n)}(x_0)$?

$newcommand{bbx}[1]{,bbox[15px,border:1px groove navy]{displaystyle{#1}},}
newcommand{braces}[1]{leftlbrace,{#1},rightrbrace}
newcommand{bracks}[1]{leftlbrack,{#1},rightrbrack}
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]{leftvert,{#1},rightvert}$

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}

For the derivation with the Newton polynomial, check this post, where the Newton polynomial is derived using a discrete calculus approach and then it is used to find the Taylor Series.

For a derivation with repeated integration, go to this site here.