% this has the integration for Simpsons rule, n must be odd
%
% [area,h,error] = simp(a,b,n)

function [area,h,error] = simp(a,b,n)

if mod(n,2) ~= 1
  disp('n must be odd, try again!')
  disp(' ')
  break
end 

% interval length:
h = (b-a)/(n-1);

for i=1:n
% n mesh points, n-1 intervals:
  x(i) = a + (i-1)*h;
end

area = 0;

for i=1:(n-1)/2
% have to be careful here.  The first step is the interval [x(1),x(3)]
% with midpoint x(2), the final step is the interval [x(n-2),x(n)] with
% midpoint x(n-1), as desired.
  j = 2*i;
  area = area + f(x(j-1)) + f(x(j))*4 + f(x(j+1));
end
area = area*h/3;

error = (f_int(b)-f_int(a))-area;