Let $Y: \mathbb{R} \to \mathbb{R}$ be three times continuously differentiable, let $x_n < x_{n+1}$ be real numbers, and let $h := x_{n+1} - x_n$. Is it true that
$$Y(x_{n+1}) = Y(x_n) + \frac{h}{2}\bigl[Y'(x_n) + Y'(x_{n+1}) \bigr]\\ - \frac{h^3}{24} \bigl[Y^{(3)}(x_n)+Y^{(3)}(x_{n+1})\bigr]+O(h^5) \quad ? \qquad (1) $$
This question arose from an answer to a previous post of mine. Formula (1) looks a lot like the Trapezoidal rule: \begin{equation} \hspace{1.5cm} \int_{x_{n}}^{x_{n+1}} f(x) \,dx = \frac{h}{2}\left[f(x_n) + f(x_{n+1}) \right] - \frac{h^3}{12} f''(\eta), \quad \eta \in [x_n,x_{n+1}]. \hspace{1.5cm} (2) \end{equation}
Indeed, if $Y$ is an anti-derivative of $f$, then we can write (2) as \begin{align*} \hspace{1.5cm} Y(x_{n+1}) = Y(x_n) + \frac{h}{2}\left[Y'(x_n) + Y'(x_{n+1}) \right] - \frac{h^3}{12} Y^{(3)}(\eta). \hspace{1.5cm} (3) \end{align*}
Comparing (1) and (3), we see that it suffices to show that \begin{equation} - \frac{h^3}{24} \bigl[Y^{(3)}(x_n)+Y^{(3)}(x_{n+1})\bigr]+O(h^5) = - \frac{h^3}{12} Y^{(3)}(\eta), \end{equation} or equivalently, for any continuous function $g$, \begin{equation} \frac{1}{2}\bigl[g(x_n)+g(x_{n+1})\bigr]+O(h^2) = g(\eta) \end{equation} for some $\eta \in [x_n,x_{n+1}]$. But I'm having trouble proving the last statement. I think Taylor series might be a useful alternative approach, but I haven't been able to make much progress with this approach...Any help would be appreciated.
You misunderstood the direction of the argument. It is the symmetry of the situation that allows to give a symmetric error term that has order gap 2 to the next error term. Using an intermediate-value argument this can then be reduced to the given error term, losing the symmetry. In the other direction the difference between the value at some midpoint and the mean value is just $O(h)$, which would still be sufficient for the claim in the previous question.
Using a double step size we know that symmetrization extracts the even and odd parts of the Taylor expansion, $$ \frac{y(x+h)-y(x-h)}{2h}=\sum_{k=0}^\infty\frac{h^{2k}}{(2k+1)!}y^{(2k+1)}(x)\\~\\ \frac{y^{(m)}(x+h)+y^{(m)}(x-h)}2=\sum_{k=0}^\infty\frac{h^{2k}}{(2k)!}y^{(2k+m)}(x) $$ Now one can use the second formula with $m=1,3,5,...$ to eliminate the derivatives on the right side in the first formula, from the lowest order on.
$$ \frac{y(x+h)-y(x-h)}{2h}-\frac{y'(x+h)+y'(x-h)}2 =-\sum_{k=1}^\infty\frac{2kh^{2k}}{(2k+1)!}y^{(2k+1)}(x) $$
$$ \frac{y(x+h)-y(x-h)}{2h}-\frac{y'(x+h)+y'(x-h)}2+h^2\frac{y'''(x+h)+y'''(x-h)}6 \\=-\sum_{k=0}^\infty\frac{h^{2k+2}}{(2k+3)(2k+1)!}y^{(2k+3)}(x) +\sum_{k=0}^\infty\frac{h^{2k+2}}{3(2k)!}y^{(2k+3)}(x) \\ =\sum_{k=1}^\infty\frac{((2k+1)(2k+3)-3)h^{2k+2}}{3(2k+3)(2k+1)!}y^{(2k+3)}(x) \\ =\sum_{k=0}^\infty\frac{(2k+6)h^{2k+4}}{3(2k+5)(2k+3)(2k+1)!}y^{(2k+5)}(x) $$
Finally, replacing $h$ back to $h/2$, $x$ to $x_{n+1/2}=x_n+\frac h2$ gives $$ \frac{y(x_{n+1})-y(x_n)}{h}-\frac{y'(x_{n+1})+y'(x_n)}2+h^2\frac{y'''(x_{n+1})+y'''(x_n)}{24} \\ =\frac{h^{4}}{120}y^{(5)}(x)+\frac{h^{6}}{5040}y^{(7)}(x)+...=O(h^4) $$ This is the claimed formula, with some extra terms.