function I=simpson(f,a,b,n)
I=2/3*midpoint(f,a,b,n)+1/3*trap(f,a,b,n);