Here is an example of Simpson's rule on (x0,y0), (x1,y1), (x2,y2) with x_0, x_0+h, x_0+2h (with y0,y1,y2). We use quadratics. restart; x1:=x0+h; x2:=x0+2*h;

LCZJI3gwRzYiIiIiSSJoR0YkRiU=

LCZJI3gwRzYiIiIiSSJoR0YkIiIj

Here are the Lagrange Polynomials: L0:=(x-x1)*(x-x2)/((x0-x1)*(x0-x2));

LCQqKCwoSSJ4RzYiIiIiSSN4MEdGJiEiIkkiaEdGJkYpRicsKEYlRidGKEYpRiohIiNGJ0YqRiwjRiciIiM=

L1:=(x-x0)*(x-x2)/((x1-x0)*(x1-x2));

LCQqKCwmSSJ4RzYiIiIiSSN4MEdGJiEiIkYnLChGJUYnRihGKUkiaEdGJiEiI0YnRitGLEYp

L2:=(x-x0)*(x-x1)/((x2-x0)*(x2-x1));

LCQqKCwmSSJ4RzYiIiIiSSN4MEdGJiEiIkYnLChGJUYnRihGKUkiaEdGJkYpRidGKyEiIyNGJyIiIw==

The error term is a function of x (below, F(x) is standing in for f'''(c(x))/3!) E:=F(x)*(x-x0)*(x-x1)*(x-x2);

KiotSSJGRzYiNiNJInhHRiUiIiIsJkYnRihJI3gwR0YlISIiRigsKEYnRihGKkYrSSJoR0YlRitGKCwoRidGKEYqRitGLSEiI0Yo

The full function is P. Differentiate twice and evaluate at x0: P:=y0*L0+y1*L1+y2*L2;

LCgqKkkjeTBHNiIiIiIsKEkieEdGJUYmSSN4MEdGJSEiIkkiaEdGJUYqRiYsKEYoRiZGKUYqRishIiNGJkYrRi0jRiYiIiMqKkkjeTFHRiVGJiwmRihGJkYpRipGJkYsRiZGK0YtRioqKkkjeTJHRiVGJkYyRiZGJ0YmRitGLUYu

P1:=simplify(int(P,x=x0..x0+2*h));

LCQqJkkiaEc2IiIiIiwoSSN5MEdGJUYmSSN5MUdGJSIiJUkjeTJHRiVGJkYmI0YmIiIk

int(E,x=x0..x0+2*h);

LUkkaW50RzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliRzYiNiQqKi1JIkZHRic2I0kieEdGJyIiIiwmRi1GLkkjeDBHRichIiJGLiwoRi1GLkYwRjFJImhHRidGMUYuLChGLUYuRjBGMUYzISIjRi4vRi07RjAsJkYwRi5GMyIiIw==

The error term comes from the MVT for integrals, so Maple does not simplify it.