restart;
# Find Newton-Cotes Quadrature Rules
n:=3;
x:=i->i/n;
L:=k->product((X-x(i)),i=0..n)/(X-x(k));
p:=add(L(k)*f(x(k))/subs(X=x(k),L(k)),k=0..n);
int(p,X=0..1);