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