newtoncotes(n) := (
  local(p,i,z0,lp,df),
  p:[],
  for i:0 thru n do p:append(p,[ [z0+i*h,y[i]] ]),
  lp : ratsimp(lagrange(p)),
  df : integrate(''lp, x, z0, z0+n*h),
  factor(ratsimp(df))
);