Consider a function, $f(x)$, which doesn't change too rapidly over the interval, or more accurately we limit the functions this method will work on to those that have no “infinities” over the interval. Examples of bad candidates for this method are the one-over functions of trigonometry with an interval of a full period (or the select part), and since we're also considering sequential derivatives also no points with an infinity in a derivative over a discontinuity (step function), on the interval, $[x_1, x_2]$. Let the $n$'th derivative of $f$, which can be written, $(d/dx)^nf=(D_x)^nf=D_x^n f$, be nonzero, implying all preceding derivatives are also nonzero, where a function must be zero everywhere to be considered a zero function (a derivative is a function with an implicit relationship to its integrated function). Then we can integrate such function, $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 the terms (adding and subtracting) 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.