%
% Composite Simpson's rule 
% In this example, we observe the accuracy of the rule 
% with respect to the number of subintervals. 

% definition of integrand
f = inline('1/x') 
a = 1;
b = 2; 
exactval = log(2); 

% loop on n = number of subintervals 
k = 0; 
for n = 2:2:100, 
   k = k+1;
   h(k) = (b-a)/n;  

   % initialization 
   sum(k) = f(a) + 4*f(a+h(k)) + f(b);

   % loop of composite Simpson's rule 
   for i = 2:2:n-2,
    sum(k) = sum(k) + 2*f(a+i*h(k)) + 4*f(a+(i+1)*h(k));
   end

   sum(k) = h(k)*sum(k)/3; 

   error(k) = abs(exactval - sum(k))/abs(exactval); 

end 

compositesimpson = sum(k)  
exactval 
relerror = error(k)

semilogy(error,'r*')
title('Composite Simpson rule') 
xlabel('n/2') 
ylabel('Relative error') 
