Consider Taylor polynomial
$$P_{n,a}(x) = \sum_{k=0}^{n} \frac{f^{(k)}(a)}{k!}(x-a)^k$$
Approximating $f(x)$ in a neighborhood of $a$. Then we have remainder function with the relationship $P_{n,a}(x) + R_{n,a}(x) = f(x)$ or $R_{n,a}(x) = f(x) - P_{n,a}(x)$.
Take $n+1$ derivatives of both sides as we get $R^{(n+1)}_{n,a}(x) = f^{(n+1)}(x)$.
At this point we're supposed to bound the derivative so when we integrate back up we back error bounds for $R_{n,a}(x)$. But does it even make sense to write
$$|f^{(n+1)}(x) | < M$$
I am a little stuck because apparently we can only bound the error between $x$ (the point we care about) and $a$ (some nearby value for which we can easily compute $f^{(k)}(a)$ values). But we already have an $x$ in the term?
I mean I don't think I can say $|f^{(n+1)}(x) | < M$ $\forall x \in (a, x)$ because that doesn't make sense. Or am I supposed to write $|f^{(n+1)}(t) | < M$ $\forall t \in (a, x)$? Like what's the usual way people write this so that when they do the $n+1$ integrations they end up at something that looks like $|R_{n,a}(x)| < \frac{M (x-a)^{n+1}}{(n+1)!}$ or whatever it's supposed to be? What's the correct range of integration? $a$ to $x$? What about when $a > x$, so $x$ to $a$? Do we have to prove both cases? Or do you do generic bounds?
I wish there were clearer information about this online in its derivation. Most websites skip all this stuff and jump right to the end result with very little explanation in between.
It is true that $R_{n,a}^{(n+1)}(x)=f^{(n+1)}(x)$. We even know $R_{n,a}^{(k)}(a)=0$ for $k=0,1,\dots,n$. Thus we actually have a well-defined ODE IVP for $R_{n,a}$, if we are given a function $f$ which is $(n+1)$ times continuously differentiable. Nonetheless this is not generally possible to solve explicitly, so it is not really useful by itself.
On the other hand, one can use this as an alternative derivation of a Lagrange-type error estimate, since it implies that the error is exactly given by:
$$R_{n,a}(x)=\int_a^x \int_a^{t_1} \dots \int_a^{t_n} f^{(n+1)}(t_{n+1}) dt_{n+1} dt_n dt_{n-1} \dots dt_1.$$
Taking absolute values and applying the triangle inequality results in the usual Lagrange-type error estimate $|R_{n,a}(x)| \leq M \frac{(x-a)^{n+1}}{(n+1)!}$ where $M$ is a bound on $|f^{(n+1)}(t)|$ for $t \in (a,x)$.
The actual Lagrange error formula says something a bit better than this estimate, namely that the error is exactly given by $f^{(n+1)}(\xi(x)) \frac{(x-a)^{n+1}}{(n+1)!}$ for some $\xi(x) \in (a,x)$. But since $\xi(x)$ generally cannot be computed, most of the time the Lagrange error formula is used to obtain the Lagrange type error estimate.