nmin=2;
nmax=100;


deff('y=f(x)','y=abs(cos(3.7*x))^2+cos(x)');

integral=zeros(nmax-nmin+1,4);

nn=3000;
h=1/nn;
xh=linspace(0,1-h,nn);
intexact=h*1/2*(sum(f(xh+h*(1/2-1/(2*sqrt(3)))))+sum(f(xh+h*(1/2+1/(2*sqrt(3))))));


for n=1:(nmax-nmin+1)
	nn=nmin+n-1;
	h=1/nn;
	xh=linspace(0,1-h,nn);
	integral(n,1)=h*sum(f(xh)); //rectangles
	integral(n,2)=integral(n,1)-h/2*f(0)+h/2*f(1); //trapezes
	integral(n,3)=1/3*integral(n,2)+2/3*h*sum(f(xh+h/2)); //Simpson
	integral(n,4)=h*1/2*(sum(f(xh+h*(1/2-1/(2*sqrt(3)))))+sum(f(xh+h*(1/2+1/(2*sqrt(3)))))); //Gauss
end;

xbasc();
x=nmin:nmax;x=log(x);
x=[x;x;x;x];
intlog=abs(integral-intexact);
intlog=log(intlog+10^(-15));
plot2d(x',intlog,leg='rectangles@trapezes@Simpson@Gauss');
