Simpson's Rule Integration

function I = simpson(a,b,N)

% This function numerically integrates a function f(x), defined in the
% function "f.m"

% It takes input a = lower bound of integration, b = upper bound of
% integration, N = number of "rectangles" we are finding the area of

sum_simp = 0; % initializing the sum of the Simpson terms
dx = (b-a)/(N); % the step between sections

for r = 1:N+1 % there are N sections (rectangles), and so there are N+1 points

    X = a + ((r-1)*dx); % each point expressed in terms of a, r and dx
    Y = f(X); % the function we are integrating

    if (mod(r,2) == 0) % essentially, "if r is even"
        coeff = 4; % setting the coefficient of all "even" terms to 4 (as per Simpson's 1/3 rule)
    else % "if r is not even"
        coeff = 2; % setting the coefficient of all "odd" terms to 2 (as per Simpson's 1/3 rule)
    end
    if (r == 1)||(r == N+1) % setting the coefficient of the first and last term to 1 (as per Simpson's 1/3 rule)
        coeff = 1;
    end

    sum_simp = sum_simp + (coeff*Y*dx*(1/3)); % adding each Simpson term to reach the final result
end
I = sum_simp; % the final output

disp(' '); disp('q1b)')
disp(strcat('The integral using Simpsons rule = ', num2str(I))); % displays the result
end