%composite Simpson's rule nparis of subintervals
function A=simprule(fh,a,b,npairs)
h=(b-a)/(2*npairs);
sumeven=0;
sumodd=0;
for i=2:2:2*npairs-2
    sumeven=sumeven+fh(a+i*h);
end
for i=1:2:2*npairs-1
    sumodd=sumodd+fh(a+i*h);
end
A=h*(fh(a)+4*sumodd+2*sumeven+fh(b))/3;