Let $L$ be the quadratic Lagrange interpolating polynomial for a function $f$ such that $f(-h)=L(-h)$, $f(0)=L(0)$, and $f(h)=L(h)$. Then it is well-known that $$f(x)=L(x)+\frac{f'''(\xi(x))}{6}x(x-h)(x+h),$$ where $\xi(x)\in(-h,h)$.
Integrating, we get $$\int_{-h}^h f(x)\,dx = \frac{h}{3}(f(-h)+4f(0)+f(h))+\int_{-h}^h\frac{f'''(\xi(x))}{6}x(x-h)(x+h)\,dx.$$
Of course, since $x(x-h)(x+h)$ changes signs on $[-h,h]$, we can't use the mean value theorem.
I have seen the following lemma used to fix this problem: $$\int_{-h}^h\frac{f'''(\xi(x))}{6}x(x-h)(x+h)\,dx=\frac{-f^{(4)}(\mu)}{24}\int_{-h}^hx^2(x-h)(x+h)\,dx$$ for some $\mu\in(-h,h)$.
Then the right hand side evaluates to $-\frac{f^{(4)}(\mu)}{90}h^5$ as desired.
What is the proof of this lemma? I imagine it's integration by parts followed by mean value theorem, but I can't seem to figure it out.