function sm = smpsn(f,rng,N)

M = N/2;
if (round(M)~=M)
    error('N must be even')
end
a = rng(1); b = rng(2);
h = (b-a)/N;
x = a + h*(1:2:(N-1));
y = a + h*(2:2:(N-2));
sm = h*(f(a)+f(b)+4*sum(f(x))+2*sum(f(y)))/3;