clear;clc;
N = 5;
f = @(x) sqrt(exp(x)-1) - x.^0.5 - x.^1.5/4 - 5*x.^2.5/96;
a = 2/3 + 1/10 + 5/336;
for i = 1:N
  n = 100*2^(i-1);
  h = 1/n;
  x = [0:h:1];
  index = 2:2:n;
  C = 1/sqrt(3);
  I(i) = h*sum(f(x(index)+C*h) + f(x(index)-C*h));
end
I = I + a;
o = (I(1:N-2)-I(2:N-1))./(I(2:N-1)-I(3:N))