function s = comp_simp(f,a,b,m)

%    A MATLAB CODE PERFORMING SIMPSON'S INTEGRATION
%
% INPUT:     f = the integrand input as a string; 
%                if the integrand is, say, f(x) = sin(x), 
%                give as first argument of comp_simp the expression
%                        inline('sin(x)')
%            a = lower limit of integration
%            b = upper limit of integration
%            m = number of intervals in which [a,b] is divided 
%                NOTE:  m must be an EVEN positive integer!
%
%            The stepsize h is equal to  (b-a)/m 
%
% OUTPUT:    The value of the definite integral

format long

h = (b-a)/m
sum_odd  = 0.0;
sum_even = 0.0;

for k = 1:(m/2)
  x = a + h * (2*k-1);
  sum_odd = sum_odd + feval(f,x);
end

for k = 1:(m/2-1)
  x = a + h * 2 * k;
  sum_even = sum_even + feval(f,x);
end

s = h * ( feval(f,a) + feval(f,b) + 4 * sum_odd + 2 * sum_even ) / 3.0;
