function adapt_quadr
a=0;     % the interval of integration [a,b]
b=pi;
e=0.00004;
%%%%%
format long
[Q, nodes]=quadr(a,b,e);
int_exact(a,b)-Q
nodes;
size(nodes)

%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
x=[a:(b-a)/100:b];
plot(x,f(x),'-r',nodes,nodes.*0,'xb');

function [S nodes]=quadr(a,b,e)
Err = abs(simps(a,b) - simps(a,(a+b)/2) - simps((a+b)/2, b));
if Err < 15*e
 S = simps(a,(a+b)/2) + simps((a+b)/2, b);
 nodes = [a, 3*a/4+b/4];
else  
 [S1, nodes1] = quadr(a,(a+b)/2,e/2); 
 [S2, nodes2] = quadr((a+b)/2,b,e/2);   
S=S1+S2;
nodes = [nodes1 nodes2];   % add the new left endpoint to check which intermediate points are introduced 

end

function S=simps(a,b);
S=(b-a)*(f(a)+4*f((a+b)/2)+f(b))/6;

function y=f(x)
%y=1/sqrt(1-x^2);  % int_{a}^{b}\frac{1}{\sqrt{1-x^2}} dx 
%y=cos(x);
y=x.*cos(x.^2);

function y=int_exact(a,b);
%y=asin(b)-asin(a);
%y=sin(b)-sin(a);
y=sin(b^2)/2-sin(a^2)/2;