For $h > 0$ consider the central difference formula
$$D(x_0) = \frac{f(x_0 + h) - f(x_0-h)}{2h} \approx f'(x_0).$$ If $f \in C^3[x_0-h, x_0+h]$, then one can show that $$|f'(x_0) - D(x_0)| \leq \frac{M}{6}h^2, \quad \text{where} \quad M = \max_{x_0 - h \leq x \leq x_0 + h} f'''(x).$$
Namely, use Taylor's formula with the Lagrange form of the remainder:
$$f(x_0 + h) = f(x_0) + h f'(x_0) + \frac{h^2}{2} f''(x_0) + \frac{h^3}{6} f'''(x_0 + \theta_1 h),$$ $$f(x_0 - h) = f(x_0) - h f'(x_0) + \frac{h^2}{2} f''(x_0) - \frac{h^3}{6} f'''(x_0 + \theta_2 h),$$ with $0 < \theta_1, \theta_2 <1$. By taking the difference of the two expressions, dividing by $2h$ and rearranging, one arrives at
$$ \frac{f(x_0 + h) - f(x_0-h)}{2h} -f'(x_0) = \frac{h^2}{12}(f'''(x_0 + \theta_1 h)+f'''(x_0 + \theta_2 h)) $$ which yields the result.
Now my question concerns another proof that I have come across which considers $$f(x_0 + h) = f(x_0) + h f'(x_0) + \frac{h^2}{2} f''(x_0) + r_3,$$ $$f(x_0 - h) = f(x_0) - h f'(x_0) + \frac{h^2}{2} f''(x_0) - R_3.$$ Then the claim is that
$$r_3 = \frac{1}{6} \int_{x_0}^{x_0 + h} (f''(x) - f''(x_0))(x_0 + h - x) dx,$$ $$R_3 = \frac{1}{6} \int_{x_0}^{x_0 - h} (f''(x) - f''(x_0))(x_0 - h - x) dx.$$
Could someone clarify why this is so and why does it yield the result?
EDIT: I tried to proceed using the formula derived by Lutzl
So we have
$$f(x_0 + h) = f(x_0) + h f'(x_0) + \frac{h^2}{2} f''(x_0) + r_3,$$ $$f(x_0 - h) = f(x_0) - h f'(x_0) + \frac{h^2}{2} f''(x_0) - R_3,$$
with
$$r_3 = \int_{x_0}^{x_0 + h} (f''(x) - f''(x_0))(x_0 + h - x) dx,$$ $$R_3 = \int_{x_0}^{x_0 - h} (f''(x) - f''(x_0))(x_0 - h - x) dx.$$
Consider the difference
\begin{align} &f(x_0 + h) - f(x_0 - h) - 2h f'(x_0) \\ &= \int_{x_0}^{x_0 + h} (f''(x) - f''(x_0))(x_0 + h - x) dx + \int_{x_0}^{x_0 - h} (f''(x) - f''(x_0))(x_0 - h - x) dx \\ &= \int_{x_0}^{x_0 + h} f'''(\xi (x))(x - x_0)(x_0 + h -x) dx + \int_{x_0}^{x_0 - h} f'''(\eta (x))(x - x_0)(x_0 - h - x) dx \\ &= \int_{x_0}^{x_0 + h} f'''(\xi (x))(x - x_0)(x_0 + h -x) dx + \int_{x_0 - h}^{x_0} f'''(\eta (x))(x - x_0)(x + h - x_0) dx \end{align}
Now divide by $2h$ and take the absolute value
\begin{align} &\left|\frac{f(x_0+h) - f(x_0 - h)}{2h} - f'(x_0)\right| \\ &\leq \frac{1}{2h} \left[ \int_{x_0}^{x_0 + h} M h (x_0 + h -x)dx + \int_{x_0 -h}^{x_0} M h (x + h -x_0)dx \right] \\ &= \frac{M}{2} \left[ \frac{h^2}{2} + \frac{h^2}{2}\right] = \frac{M}{2} h^2 \end{align}
So the esitmate is $\frac{M}{2} h^2$ unlike the initial $\frac{M}{6} h^2$. Does this seem right?
From partial integration you know that \begin{align} f(x+h)&=f(x)+f'(x)h+\int_x^{x+h}(x+h-s)f''(s)\,ds \end{align} Now $$ \int_x^{x+h}(x+h-s)\,ds=\frac12h^2 $$ so that adding a zero we get \begin{align} f(x+h)&=f(x)+f'(x)h+\frac12f''(x)h^2+\int_x^{x+h}(x+h-s)(f''(s)-f''(x))\,ds. \end{align} It is not clear where the different constants in the given formula come from. For functions with $f''(x)=0$ they give a visibly wrong result.
Now if you know that $|f''(x)-f''(y)|\le M|x-y|$, then (using $s=x+th$) $$ \left|\int_x^{x+h}(x+h-s)(f''(s)-f''(x))\,ds\right|\le M\int_x^{x+h}(x+h-s)(s-x)\,ds =Mh^3\int_0^1(1-t)t\,dt=\frac16Mh^3. $$ If you insert this into the error estimate formula, you get $\frac16Mh^2$ as bound.