function Inter3 = SimpsonRule(f,a,b,n)
% This function integrates using the Simpson Rule where:
% f = function handle
% a = lower limit (start)
% b = upper limit (end)
% n = number of sections

% This first section is implemented incase of multiple values of n
numbofsections = n;  
numbofn = size(n);
column = numbofn(1,2);
Inter3(1,column) = NaN;
for column = 1:column
n = numbofsections(1,column);
h = (b-a)./n;
x=linspace(a,b,n+1);
y_val(1,n+1) = NaN;

% This for loop calculates values of y from f(x)
for i=1:n+1
y_val(1,i) = feval(f,x(i));
end

% Indexing and retrieving odd numbers from y-vals
oddindex = @(y_val) y_val(1+1:2:end-1);
yodd = oddindex(y_val);

% Indexing and retrieving even numbers from y-vals
evenindex = @(y_val) y_val(1+2:2:end-1);
yeven = evenindex(y_val);

% This line is the formula for calculating the estimated area.
Inter3(1,column) = (h/3).*((y_val(1)+y_val(end))+(4.*sum(yodd))...
    +(2.*sum(yeven)));

end
end