Consider a function, $f(x)$, which doesn't change too rapidly over the interval, or we limit the functions this method will work on to those which have no “infinities” on the interval. For example, the one-over functions of trigonometry with an interval of a full period (or the select part) have infinities, because they hold at least the first-order singularity $\lim_{x\rightarrow 0}\frac{1}{x}$.
Let the $n$'th derivative of $f$ exist (that is, be non-zero), written, $(d/dx)^nf=(D_x)^nf=D_x^n f$, implying all preceding derivatives are also nonzero, where a function must be zero everywhere to be considered a zero function. For example, the exponential and trigonemetric functions all have non-zero derivatives, whereas a polynomial only has non-zero derivatives up to the highest order term. Then we can integrate such a Taylor polynomial, $F=D_x^n f$, anywhere inside the interval, say $[x_1, x]$, where $x\lt x_2$, and consider the results.
$$ I_1 = \int_{x_1}^{x} F dx = \int_{x_1}^{x} D_x^{n} f dx = \left. D_x^{n-1} f \right|_{x_1}^x = D_x^{n-1} f(x) - D_x^{n-1} f(x_1) $$Since the last term is constant and recurring, I'll label it, $c_m=-\left. D_x^{n-m}f(x) \right|_{x_1}$. Each of the $c$ terms are the result of one iteration of integrating the function $F$, the $n$th derivative of $f$, and evaluating such ($m$th) integral at the two points of the finite integral limits, $x_1$ on the left and $x\lt x_2$ on the right, with the first being subtracted from the second.
$$ I_2 = \int I_1 dx = \int_{x_1}^{x} D_x^{n-1} f dx + \int_{x_1}^x c_1 dx = D_x^{n-2} f(x) + c_2 + c_1 (x-x_1) $$ $$ I_3 = \int I_2 dx = \int_{x_1}^{x} D_x^{n-2} f dx + \int_{x_1}^x [ c_1 ( x - x_1 ) + c_2 ] dx $$ $$ = D_x^{n-3}f(x) + c_3 + c_2(x-x_1) + \frac{c_1}{2}( x^2 - 2 x_1 x - x_1^2 + 2 x_1^2 ) $$The coefficient of $c_1$ remains as the four terms of the integral, but can be compactified as $(x-x_1)^2/2$. Then continuing the integrations, the $m=n$'th step will have the following terms: $$ I_n = \int I_{n-1} dx = f(x) +c_n + \frac{c_{n-1}}{1!} (x-x_1) $$ $$ + \frac{c_{n-2}}{2!} (x-x_1)^2 +... + \frac{c_1}{(n-1)!} (x-x_1)^{n-1} $$ Where $f(x)=D^0_x f(x)$, and with some algebraic rearangement (adding and subtracting) the terms in $I_n$ yields the following, with $f(x)$ on the left, and the remainder on the right, with the $c$ terms substituted out, to obtain the final form of the Taylor series.
$$ f(x) = f(x_1) + D_x f(x_1) (x-x_1) + D_x^2 f(x_1) \frac{1}{2!}(x-x_1)^2 $$ $$ + D_x^3 f(x_1) \frac{1}{3!}(x-x_1)^3 +...+-I_n $$The remainder term, $I_n=\int_{x_1}^x \cdots \int_{x_1}^x (D_x)^n f (x) dx^n$, has an integrand whose dependence on the number of terms, $n$, is small for some functions, constant for the exponential-linear ones, and alternating for sinusoidal and others, but we can further apply analysis on the integrated term, by bounding the first integral. We can say that $F(x)$, is less than, or equal to, the magnitude of the interval, $(x-x_1)$, times the maximum of $F=D_x^n f$ on the interval, which exists since we are only talking about functions, $f$, which are well behaved on that interval. Further, we also know that $F$ has a minimum, out of all the values it may take, $F(x)$, on the interval. So that the first integral is in the ballpark of the product of the interval, $x$, times the average of $F$, $F(x_{\text{avg}})$. And so we jump ahead on calculating the subsequent steps, with an error value depending invariably upon $n$, as it diminishes the remainder term by $n!$, on the denominator, which is arbitrarily large.
$$ I_n = \frac{(x-x_1)^{n}}{n!} D_x^n f(x_{\text{avg}}) $$Immediately we can see that the approximation for a function, $f$, contains a geometric series component that tends to converge it for intervals less than unity, $(x-x_1)\lt 1$, and since the other coefficient is monotonically decreasing with $n$, as well, $n^{-1} \gt (n+1)^{-1} , \forall\ n\gt 0$, the series converges, and the amount of deviation from the true function is controlled.