function I= Simpson( a, b, n, f )
% Calcul integrala prin metoda Simpson
% f functie, a, b capetele intervalului, n precizia

    m =(b-a) / (2*n);
    s1 = 0;
    s2 = 0;
    for i = 1 : n
        s1 = s1 + f(a+(2*i-1)*m); %calculul f1
    end
    for i = 1 : n-1
        s2 = s2 + f(a+2*i*m); %calculul f2
    end
    I = m*(f(a) + f(b) + 4*s1 + 2*s2)/3;

end

% P.S. Desi m-am inspirat, am inteles ce e cu formula :D