function I = simpsons(a, b, f, N);

if mod(N, 2)
   fprintf('Error: subinterval number (N=%d) must be even\n', N);
end

dx = (b - a) / N;
S = 0;

for i = 1 : (N - 1)
       x = a + i * dx;
       S = S + (1 + mod(i, 2)) * f(x);  % Weight sequence 2, 1, 2, 1, ...
end

I = (f(a) + 2 * S + f(b)) * dx / 3;     % ...multiplied by two here
    
       