Adaptive Simpson's Rule with true error for sin(x)

function [outI, outX] = adaptive_two_simpson_f_error (func, a, b, tol)
m = (a+b)/2;
## c = 1/(1/16-1);
## k = 1/15;
## f = sin(x) for reference
I = cos(a)-cos(b);

[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(I - (outI_1 + outI_2)) < tol
  outI = outI_1 + outI_2;
  outX = unique([outX_1, outX_2]);

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

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

endif

endfunction