Simpson Three Eight

function IS = simpson_three_eight(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, 3) == 0
            sum = sum + 2 * f(a + i * h);
        else
            sum = sum + 3 * f(a + i * h);
        end

    end

    IS = (3 * sum * h) / 8;
end