CCCPaste

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