function comp_simpson
n=2*28; % the number of nodes: MAKE SURE n IS EVEN! 
a=0; b=pi; % the interval [a,b]
h=(b-a)/n; % the step size
nodes=a:h:b; % the nodes
%%%%% THE COMPOSITE SIMPSON RULE %%%%%
%%%%% TO CALCULATE \int_{a}^{b} f(x)dx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S=f(a);
for i=0:1:fix(n/2)-2
  S=S+4*f(a+i*2*h+h)+2*f(a+i*2*h+2*h);
end;
S=S+4*f(b-h)+f(b);
S=(h/3)*S;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compare to the exact integral 
%S+(cos(b)-cos(a))
S-sin(b^2)/2+sin(a^2)/2  



function y=f(x);
%y=sin(x);
y=x*cos(x^2);