Modified Composite Quadrature formula I'm trying to produce a modified composite trapezoid quadrature formula for $\int_a^b f(x)\ dx$ based on the Euler-MacLaurin formula. 
I know that it should involve $f'(a), f'(b), f'''(a)$ and $f'''(b)$, but I really don't know how to put it together. Also, I think it has an error on the order of $h^6$, but I'm having a hard time seeing the math.
 A: My favorite approach to these sorts of formulas is to note that integration is linear, and so on polynomials, this amounts to a little linear algebra. There are six variables: $f^{(n)}(a):=a_n$ and $f^{(n)}(b):=b_n$ for $n\in\{0,1,3\}$. Thus we can expect to be able to find the first six terms of the power series, and so we should be able to integrate any fifth-order polynomial exactly. So now let's apply our formula to an arbitrary such polynomial:
$$f(x)=\sum_{n=0}^5c_nx^n$$
For convenience, we'll set $a=0$, $b=1$ in the integral to simplify our equations. The integral of $f$ is $\int_0^1f(x)\,dx=\sum_{n=0}^5\frac{c_n}{n+1}$, and the derivatives of $f$ are $a_n=c_nn!$ for $n\in\{0,1,3\}$ and
$$b_0=\sum_{n=0}^5c_n\qquad b_1=\sum_{n=1}^5nc_n\qquad b_3=6c_3+24c_4+60c_5.$$
This is a linear system of six equations relating $c_n$ to $a_n,b_n$. Solving them yields:
$$c_0=a_0\qquad c_1=a_1\qquad c_2=-\frac52a_0-\frac74a_1-\frac{a_3}{16}+\frac52b_0-\frac34b_1+\frac{b_3}{48}$$
$$c_3=\frac{a_3}6\qquad c_4=\frac52a_0+\frac54a_1-\frac7{48}a_3+\frac54b_1-\frac52b_0-\frac{b_3}{16}$$
$$c_5=-a_0-\frac{a_1}{2}+\frac{a_3}{24}+b_0-\frac{b_1}{2}+\frac{b_3}{24}$$
and plugging them into the formula for the integral of $f$ (and substituting back the definitions of $a_n,b_n$) gives
$$\int_0^1f(x)\,dx=\frac{f(1)+f(0)}{2}-\frac{f'(1)-f'(0)}{12}+\frac{f'''(1)-f'''(0)}{720}$$
which is our final formula. If we are taking a general integral over $[a,b]$ instead of $[0,1]$, we can "stretch" the formula by scaling by $b-a$ horizontally. This will increase the integral by $b-a$, divide the first derivative by $b-a$, and divide the third derivative by $(b-a)^3$, so:
$$\int_a^bf(x)\,dx=\frac{f(b)+f(a)}{2}(b-a)-\frac{f'(b)-f'(a)}{12}(b-a)^2+\frac{f'''(b)-f'''(a)}{720}(b-a)^4$$
and the error term is $\int_a^b\sum_{n=6}^\infty c_nx^n\,dx=O((b-a)^6)$.

It is asked how this all relates to the Euler-Maclaurin formula, which is in this context
$$\sum_{n=0}^1f(n)=\int_0^1f(x)\,dx+\frac{f(1)+f(0)}2+\sum_{k=1}^\infty\frac{B_{2k}}{(2k)!}\left(f^{(2k-1)}(1)-f^{(2k-1)}(0)\right).$$
Truncating the sum at $k=2$ yields:
$$\int_0^1f(x)\,dx=\frac{f(1)+f(0)}2-\frac{B_2}{2!}(f'(1)-f'(0))-\frac{B_4}{4!}(f'''(1)-f'''(0)),$$
which coincides with our earlier formula because $\dfrac{B_2}{2!}=\dfrac{1/6}{2}=\dfrac1{12}$ and $\dfrac{B_4}{4!}=\dfrac{-1/30}{24}=\dfrac1{720}$.
