function main
    global fcount
	format compact
    format long
    a = -1
    b =  1
	f = @my_integrand

	fcount = 0;
	reltol = 1.0e-10;
    m0 = 8;
	mmax = 2^10;
	approx1 = MyCompositeSimpson(f,a,b,m0)
	good_enuf = false;
	m = m0;

	while ~good_enuf
	    m = 2*m
		if m>mmax
			error('Unable to get desired accuracy');
		end
		approx2 = MyCompositeSimpson(f,a,b,m)
	    error_estimate = approx1 - approx2
		if abs(error_estimate/approx2) < reltol
			good_enuf = true;
	    end
		approx1 = approx2;
	end
	answer = approx2
	fprintf('required %i function evaluations\n',fcount);
end

function I=MyCompositeSimpson(f,a,b,m)
   h = (b-a)/(2*m);
   x = a:h:b;
   factors = ones(2*m+1,1);
   for i=2:2:2*m
       factors(i) = 4;
   end
   for i=3:2:2*m-1
       factors(i) = 2;
   end
   factors';
   fx = f(x);
   I = h/3 * (fx*factors);
end

function y = my_integrand(x)
	global fcount
	y = runge(x-0.3) + 0.2*runge(5*(x-0.6));
	fcount = fcount + numel(x);
end

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