Show that $\left|\int_{-h}^hf(x)\mathrm dx-\frac{h}3(f(-h)+f(h)+4f(0))\right|\le\frac{h^5}{90}\|f^{(4)}\|_\infty$ 
Suppose that $f\in C^4([-h,h],\Bbb R)$ with $h>0$. Show that $$\left|\int_{-h}^hf(x)\mathrm dx-\frac{h}3(f(-h)+f(h)+4f(0))\right|\le\frac{h^5}{90}\|f^{(4)}\|_\infty$$
HINT: integrate by parts and after use the mean value theorem for integrals.

I tried different approaches with no results, by example I integrate by parts, prior change of variable, using the Bernoulli polynomials and tried to setup something using two different versions of the MVT for integrals.
I tried other more direct approaches without the Bernoulli polynomials with no result. The two versions of the MVT for integrals that I know are these:

MVT for integrals (version I). Let $f,g\in C(I,\Bbb R)$ with $g\ge 0$, then
$$\int_\alpha^\beta f(x)g(x)\mathrm dx=f(\xi)\int_\alpha^\beta g(x)\mathrm dx,\quad \xi\in[\alpha,\beta]$$
MVT for integrals (version II). Let $f\in C(I,\Bbb R)$ and $g\in C^1(I,\Bbb R)$ monotone, then
$$\int_\alpha^\beta f(x)g(x)\mathrm dx=g(\alpha)\int_\alpha^\xi f(x)\mathrm dx+g(\beta)\int_\xi^\beta f(x)\mathrm dx,\quad \xi\in[\alpha,\beta]$$

Some help will be appreciated, thank you.
 A: As commented, there are several approaches. However, if you want to follow the hint in brute-force fashion, then integrate by parts four times retaining constants of integration. Solve for the constants that give the desired quadrature rule and eliminate terms where $f'$, $f''$ and $f'''$ appear.
Clearly, this is not an elegant approach here, but in other contexts the technique is useful.
For the first, step take
$$\int_{-h}^h f(x) \, dx = \int_{-h}^0 f(x) \, dx + \int_{0}^h f(x) \, dx =  \int_{0}^h f(x-h) \, dx + \int_{0}^h f(x) \, dx. $$
Perform integration by parts using $dv = dx \implies v = x + \text{const.}$ to obtain
$$\int_{-h}^h f(x) \, dx \\ = \left. (x + A_1)f(x)\right|_{0}^h + \int_{0}^h (x + A_1) f'(x-h)\,dx + \left. (x + B_1)f(x)\right|_{0}^h - \int_{0}^h (x + B_1) f'(x)\,dx  \\ = (h + B_1) f(h) + (h + A_1 - B_1) f(0) - A_1 f(-h) \\ -  \int_{0}^h (x + A_1) f'(x-h)\,dx -\int_{0}^h (x + B_1) f'(x)\,dx$$
Put $A_1 = -h/3$ and $B_1 = -2h/3$ to obtain
$$\int_{-h}^h f(x) \, dx  = \frac{h}{3} \left(f(-h) + 4 f(0) + f(h) \right) - E,$$
where
$$E =  \int_{0}^h (x + A_1) f'(x-h)\,dx  + \int_{0}^h (x + B_1) f'(x)\,dx.$$
Now integrate by parts three more times to obtain 
$$\begin{align}E &= \left.\left(\frac{1}{2} x^2 + A_1 x + A_2\right)f'(x-h) \right|_0^h \\ &+ \left.\left(\frac{1}{2} x^2 + B_1 x + B_2\right)f'(x) \right|_0^h \\ &- \left.\left(\frac{1}{6} x^3 + \frac{1}{2}A_1 x^2 + A_2x + A_3\right)f''(x-h) \right|_0^h \\ &- \left.\left(\frac{1}{6} x^3 + \frac{1}{2}B_1 x^2 + B_2x + B_3\right)f''(x) \right|_0^h \\ &+ \left.\left(\frac{1}{24} x^4 + \frac{1}{6}A_1 x^3 + \frac{1}{2}A_2x^2 + A_3x + A_4\right)f'''(x-h) \right|_0^h \\ &+ \left.\left(\frac{1}{24} x^4 + \frac{1}{6}B_1 x^3 + \frac{1}{2}B_2x^2 + B_3x + B_4\right)f'''(x)\right|_0^h  \\ &- \int_0^h  \left(\frac{1}{24} x^4 + \frac{1}{6}A_1 x^3 + \frac{1}{2}A_2x^2 + A_3x + A_4\right)f^{(4)}(x-h) \, dx \\ &- \int_0^h \left(\frac{1}{24} x^4 + \frac{1}{6}B_1 x^3 + \frac{1}{2}B_2x^2 + B_3x + B_4\right) f^{(4)}(x) \, dx \end{align}$$
Solve for $A_2, A_3, A_4, B_2, B_3, B_4$ to eliminate the boundary terms and apply the bound $\|f^{(4)} \|_\infty$ inside the integrals to obtain the desired result
$$|E| \leqslant \frac{h^5}{90} \|f^{(4)} \|_\infty.$$
A: Just for the record I will add the "ugly" brute-force solution. Integrating by parts recursively we get
$$\int_{-h}^hf(x)\mathrm dx=\int_0^h f(x)\mathrm dx+\int_{-h}^0f(x)\mathrm dx=\\=p_1(x)f(x)|_{-h}^0+q_1(x)f(x)|_0^h-\int_{-h}^0p_1(x)f'(x)\mathrm dx-\int_0^h q_1(x)f'(x)\mathrm dx=\\=\sum_{k=1}^4(-1)^{k+1}\Bigg(p_k(x)f^{(k-1)}(x)\Big|_{-h}^0+q_k(x)f^{(k-1)}(x)\Big|_{0}^h\Bigg)+\int_{-h}^0 p_4(x)f^{(4)}(x)\mathrm dx+\int_0^h q_4(x)f^{(4)}(x)\mathrm dx$$
where the $p_k$ and $q_k$ are polynomials of degree $k$ such that $[p'_k]=p_{k-1}$ and $[q_k]'=q_{k-1}$. We must set the coefficients of these polynomials to fit the desired result, that is, we need that
$$p_1(x)f(x)|_{-h}^0+q_1(x)f(x)|_0^h=\frac{h}3\Big(f(-h)+4f(0)+f(h)\Big)$$
and that the other terms, but the last integral, will be zero. Hence setting $p_1(x)=x+\frac{2h}3$ and $q_1(x)=x-\frac{2h}3$ above we get the RHS. Now we set
$$p_2(x)=\frac{x^2}2+h\frac{2x}3+h^2\frac16,\quad q_2(x)=\frac{x^2}2-h\frac{2x}3+h^2\frac16$$
$$p_3(x)=\frac{x^3}6+h\frac{x^2}3+h^2\frac{x}6,\quad q_3(x)=\frac{x^3}6-h\frac{x^2}3+h^2\frac{x}6$$
$$p_4(x)=\frac{x^4}{24}+h\frac{x^3}9+h^2\frac{x^2}{12}-\frac{h^4}{72},\quad q_4(x)=\frac{x^4}{24}-h\frac{x^3}9+h^2\frac{x^2}{12}-\frac{h^4}{72}$$
Under these conditions we have that
$$p_k(x)f^{(k-1)}(x)|_{-h}^0+ q_k(x)f^{(k-1)}(x)|_0^h=0,\quad k\in\{2,3,4\}$$
Then we have the identity
$$\int_{-h}^hf(x)\mathrm dx=\frac{h}3\big(f(-h)+4f(0)+f(h)\big)+\int_{-h}^0 f^{(4)}(x)p_4(x)\mathrm dx+\int_0^h f^{(4)}(x)q_4(x)\mathrm dx$$
Hence
$$\begin{align}\left|\int_{-h}^0 f^{(4)}(x)p_4(x)\mathrm dx+\int_0^h f^{(4)}(x)q_4(x)\mathrm dx\right|&\le\int_{-h}^0|f^{(4)}(x)p_4(x)|\mathrm dx+\int_0^h|f^{(4)}(x)q_4(x)|\mathrm dx\\&\le\|f^{(4)}\|_\infty\left(\int_{-h}^0|p_4(x)|\mathrm dx+\int_0^h|q_4(x)|\mathrm dx\right)\\&\le2\|f^{(4)}\|_\infty\int_0^h\left|\frac{x^4}{24}-h\frac{x^3}9+h^2\frac{x^2}{12}-\frac{h^4}{72}\right|\mathrm dx\\&\le\|f^{(4)}\|_\infty\frac{h^5}{90}\end{align}$$
where we used the fact that $p_4(x)=q_4(-x)$ and that $q_4(x)\le0$ in $[0,h]$.$\Box$
