function adaptive_quad()
	global xvals fvals fcount
	format long;
	fcount = 0;
	a = -1
	b = 1
	reltol = 1.e-10
	integral = quad(@my_integrand,a,b,[0,reltol])
	fprintf('Required %i function evaluations\n',fcount);
	hold off;
	figure(1);
	plot(xvals,fvals,'ro'); xlabel('x'); ylabel('f(x)');
	figure(2);
	plot(xvals,1:fcount,'bo'); xlabel('x'); ylabel('f-call order');
end

function y = my_integrand(x)
	global xvals fvals fcount

	y = runge(x-0.3) + 0.2*runge(5*(x-0.6));
	fcount = fcount + 1;
	xvals(fcount) = x;
	fvals(fcount) = y;
end

function y = runge(x)
	y = 1/(1+12*x^2);
end