function J=NewtonCotes2(m,a,b,flag,N)
J=0;
for i=1:N
    J=J+NewtonCotes('f',m,a+(b-a)*(i-1)/N,a+(b-a)*i/N,flag);
end