function result=Simpson(a,b,n,F)

% a = left boundary point of integration interval
% b = right boundary point of integration interval
% n = number of intervals
% F = function handle

h=(b-a)/n;

result=0.;

% the boundary points
result=result+F(a)+F(b);
plot(a,F(a),'*')
plot(b,F(b),'*')
plot(a,0,'+')
plot(b,0,'+')

% the interior boundary points of the intervals
for i=1:n-1
   xi=a+h*i;
   result=result+2.*F(xi); 
   plot(xi,F(xi),'*')
   plot(xi,0,'+')
end

% the midpoints of the intervals
for i=1:n
   xi=a+h/2+(i-1)*h; 
   result=result+4.*F(xi);
   plot(xi,F(xi),'*')
   plot(xi,0,'+')
end

% the multiplication factor
result=h/6*result;
