I was using the Maple software for the computation of the numerical integration error. I'm obtaining two strange properties which I'm unable to explain. First for items $(72)$ and $(73)$ below while integrating those natural polynomials as displayed from $0$ to $3\cdot h$ and $5\cdot h$ resp. I'm getting the exact coefficient of error in the Simpson's rule and the $6-$Newton-Cotes rule, together with the right order of $h$: $5$ and $7$ resp. Is this a mere coincidence that these two coefficients are the same as they should be ? How could I explain this phenomena ? Next bizarre thing is that this doesn't work for $n$ even, I'm getting just zero: $0$. How should I fix the polynomial in $(74)$ with the upper bound to obtain $-\frac{8}{945}h^7$ as given in the second snippet below in the same manner as for $(72)$ and $(73)$ ??
This question is addressed in the last two sequences of the cited text. In your last example the number of points is odd, the product of the linear factors is an odd function about the middle point of the interval, so the integral is zero.
More in detail, the integral itself is invariant under direction change, so for the most simple case of 3 points and the model case $a,a+h,a+2h=b$ one gets $$ \int_{-h}^hf(x)dx=\int_{-h}^h [p_3(x)+r(x)(x+h)x(x-h)]dx \\ =\int_{-h}^h p_3(x)dx +\int_0^h [r(x)-r(-x)](x+h)x(x-h)dx $$ with $r(x)=f[x,-h,0,h]$ an iterated divided difference, so that $r(x)-r(-x)=2xf[x,-x,-h,0,h]=2x\tilde r(x)$ with $\tilde r(x)=f[x,-x,-h,0,h]=\frac1{4!}f^{(4)}(\xi_x)$, $\xi_x\in(-h,h)$.
Inserting this then gives the remainder term as $$ 2\int_0^h \tilde r(x)(x+h)x^2(x-h)dx=2\tilde r(c)\int_0^h(x^4-h^2x^2)dx =\frac{2f^{(4)}(\xi_c)}{4!}(\tfrac15-\tfrac13)h^5 $$
So in summary, for odd numbers of points the middle factor gets squared in the computation of the coefficient of the error term.