source("adapt.m");

global ncalls; ncalls=0;

function y=f(x)
	global ncalls; ncalls++; y=1/(1+x*x);
endfunction

a=-20; b=20; eps=0.0001; acc=0.0001;

disp("easy integral:"); eps, acc, format long;
exact=atan(20)-atan(-20); exact

disp("---------------------"); disp("adapt24:");
[J,err]=adapt("f",a,b,eps,acc);
J, err, ncalls, actual_err=J-exact; actual_err

disp("---------------------"); disp("built-in Octave function quad:");
[J,ier,nfun,err]=quad("f",a,b,[acc,eps]);
J, err, nfun, actual_err=J-exact; actual_err

disp("\n\ndifficult integral:");

function y=badf(x)
	global ncalls; ncalls++; y=log(x)/sqrt(x);
endfunction

max_recursion_depth=40000; ncalls=0;
a=0; b=1; eps=0.001; acc=0.001;
eps, acc, format long;
exact=-4; exact

disp("---------------------"); disp("adapt24:");
[J,err]=adapt("badf",a,b,eps,acc);
J, err, ncalls, actual_err=J-exact; actual_err

disp("---------------------");
disp("built-in Octave function quad:");
[J,ier,nfun,err]=quad("badf",a,b,[acc,eps]);
J, err, nfun, actual_err=J-exact; actual_err
