function s = simprl(f,a,b,m)
% s = simprl(f,a,b,m)
% Quadrature using Simpson`s rule.
% f is the name of the function, input.
% a is the left endpoint, input.
% b is the right endpoint, input.
% m is the number of subintervals, input.
% s is the Simpson rule sum, output.
h = (b - a)/(2*m);
s1 = 0;
s2 = 0;
for k=1:m,
x = a + h*(2*k-1);
s1 = s1 + feval(f,x);
end
for k=1:(m-1),
x = a + h*2*k;
s2 = s2 + feval(f,x);
end
s = h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;