simpson comp

function int = simpsoncomp(f, m, a, b)
    h=(b-a)/(2*m);
    for i=1:(2*m)
        y(i)=f(a+h*i);
    end
    int = (h/3)*(f(a)+f(b)+4*sum(y(1:2:m))+2*sum(y(2:2:m-1)));
end