Simpson One Third

function IS = simpson_one_third(a, b, n)

    f = @(x) sqrt(1 - 0.162 * (sin(x) ^ 2));
    h = (b - a) / n;

    sum = f(a) + f(b);

    for i = 1 : n - 1

        if mod(i, 2) == 0
            sum = sum + 2 * f(a + i * h);
        else
            sum = sum + 4 * f(a + i * h);
        end

    end

    IS = (sum * h) / 3;
end