Composite Simpson's Rule Matlab

function I = CSR(f,a,b,n)
h=(b-a)/n;
XI=f(a)+f(b);
X_odd=0;
X_even=0;
for i=1:n-1
    if mod(i,2)==0
        X_even=X_even+f(a+i*h);
    else
        X_odd=X_odd+f(a+i*h);
    end
end
I=(XI+2*X_even+4*X_odd)*h/3;
end