Adaptive Simpson's Rule

function [outI, outX] = adaptive_two_simpson (func, a, b, tol)
m = (a+b)/2;
##c = 1/(1/16-1);
k = 1/15;

[outI_0, outX_0] = composite_two_simpson(func, a, b, (b-a)/2 );
[outI_1, outX_1] = composite_two_simpson(func, a, m, (m-a)/2 );
[outI_2, outX_2] = composite_two_simpson(func, m, b, (b-m)/2 );

if abs(k*(outI_0 - (outI_1 + outI_2))) < tol
  outI = outI_1 + outI_2;
  outX = unique([outX_1, outX_2]);

else
  [outI_3, outX_3] = adaptive_two_simpson(func, a, m, tol/2);
  [outI_4, outX_4] = adaptive_two_simpson(func, m, b, tol/2);

  outI = outI_3 + outI_4;
  outX = unique([outX_3, outX_4]);

endif

endfunction