%% Simpson's Rule - Simplest Implementation

% This code requires that the function to be integrated is defined in a
% separate file (fun.m)

%% Initialization  

clc
a = 0; % the end points of the interval of integration
b = 3;

n = 6; % MUST BE EVEN

%% Test for n even or odd
    if n/2>floor(n/2) 
        
        disp('Your n is odd! Simpson Rule requires an even n.')
        disp('Please rerun the code with an even integer n.')
    else

        dx = (b-a)/n;

        %% Numerical Integration (Simpson's Rule)

        xi = a:dx:b;

        z = fun(xi(1));

        k = n/2;

        for j=1:k

            z = z + 4*fun(xi(2*j))+2*fun(xi(2*j+1));

        end
        z = z - fun(xi(n+1));

        z = z*(b-a)/(3*n);

        disp('The approximate value of the integral (using Simpson Rule) is ')
        disp(z)

    end