function S = compsimp(a,b,N)
% Composite Simpson's Rule for function f(x) 
% on [a,b] using 2N steps
f = @(x) exp(x);
h = (b-a)/(2*N);
i = 0:N-1;
xi = a+2*i*h;
xi1 = a+2*(i+0.5)*h;
xi2 = a+2*(i+1)*h;
S = (h/3)*sum(f(xi)+4*f(xi1)+f(xi2));
end

