function I = Simpson13(f, a, b, N)

h = abs(b-a)/N;

s_par = 0;
s_impar = 0;

for k=1:N-1
    if mod(k,2) == 1
        s_impar = s_impar + f(a+h*k);
    else
        s_par = s_par + f(a+h*k);        
    end;
end;

I = (b-a)*(f(a) + 4*s_impar + 2*s_par + f(b))/(3*N);

end