% --------------------
% A general program for performing numerical
% integration using Simpson's 1/3 method
%
% (1) The function in the integrand, f(x), is defined
%     by the first line (the inline statement)
% (2) The interval for integration is [a,b]
% (3) ndeltax is the number of subintervals that 
%     the full interval is divided into; it must 
%     be an even number.
%  Written by HPH
% ---------------------
f = inline('sin(x)','x');
a = 0;
b = 1;
ndeltax = 12;
ngridpoint = ndeltax+1;
dx = (b-a)/ndeltax;
%
% create weights 
%
w(1) = 1;
w(ngridpoint) = 1;
for k = 2:ngridpoint-1
    if mod(k,2) == 0
        w(k) = 4;
    else
        w(k) = 2;
    end
end
%
% define grid points x(i)
%
for k = 1:ngridpoint
    x(k) = (k-1)*dx;
end
%
% the main integration
%
integ = 0;
for k = 1:ngridpoint
    integ = integ+f(x(k))*w(k);
end
%
% multiply the final integral by (delta-x)/3
% (the common factor for Simpson's 1/3 method),
% then output the final answer
%
integ = integ*(dx/3);
fprintf('result of integration is %8.5f \n',integ)