function IW = weddles(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 md = mod(i, 6); func_val = f(a + i * h); if md == 1 || md == 5 sum = sum + 5 * func_val; elseif md == 2 || md == 4 sum = sum + func_val; elseif md == 3 sum = sum + 6 * func_val; else sum = sum + 2 * func_val; end end IW = (3 * sum * h) / 10; end