ClosedNewtonCotesRule:=n->factor(int(interp([seq(i*h,i=0..n)],
                          [seq(f(i*h),i=0..n)],z),z=0..n*h)):