function [I,N] = simpson_composto(fun, a, b, tol, hmin)

N = 1;
S0 = feval(fun, a) + feval(fun,b);
h = (b - a)/(2 * N);
Simpar = feval(fun, a+h);
Spar = 0;
I = h * (S0 + 4 * Simpar + 2 * Spar)/3;

flag = 1;

while (flag)
    Iold = I;
    N = N + N;
    h = 0.5 * h;
    Spar = Spar + Simpar;
    Simpar = 0;
    for j = 1:N
        Simpar = Simpar + feval(fun, a + (2*j-1)*h);
    end
    I = h * (S0 + 4 * Simpar + 2 * Spar)/3;
    if abs(I - Iold) < tol | h < hmin
      flag = 0;
    end
end
I = (16 * I - Iold)/15;
