function [J,err] = adapt(f,a,b,eps,acc) 
	f2=feval(f,a+(b-a)/3); f3=feval(f,a+2*(b-a)/3);
	[J,err]=adapt24(f,a,b,f2,f3,eps,acc);
endfunction

function [J,err]=adapt24(f,a,b,f2,f3,eps,acc)
	h=b-a; f1=feval(f,a+h/6); f4=feval(f,a+5*h/6);
	j4=(f1+f4)*h/3+(f2+f3)*h/6; j2=(f1+f2+f3+f4)*h/4;
	tol=acc+eps*abs(j4); err=abs(j4-j2)/2;

	if (err<tol)
		J=j4;
	else
		[i1,err1]=adapt24(f,a,(a+b)/2,f1,f2,eps,acc/sqrt(2));
		[i2,err2]=adapt24(f,(a+b)/2,b,f3,f4,eps,acc/sqrt(2));
		J=i1+i2; err=sqrt(err1*err1+err2*err2);
	endif
endfunction
