function A=compsimp(f,a,b,m)
%Approximates the integral of f over [a,b] using 2m+1 points.

h=(b-a)/(2*m);
x=linspace(a,b,2*m+1);
y=f(x);
A1=4*sum(y(2:2:(2*m)));
A2=2*sum(y(3:2:(2*m-1)));
A=(h/3)*(y(1)+y(2*m+1)+A1+A2);