function approx = simpson(fctn, a, b, N)
%  Compute the integral of the function given by the function handle
%  fctn between the limits a and b using Simpson's Method with N 
%  subintervals.
%  N must be an even integer and b must be greater than a.
%

h = (b-a) / N;             % compute the step size for N subintervals
sum = 0;                   % intialize sum for storing results

% compute the term for each subinterval and add to sum
for n_count = 1 : N+1
    
    % get n from n_count
    n = n_count - 1;

    % compute the next x value
    x_n = a + n*h;
    
	% compute the next simpson’s value
    sum = sum + coefficient(n, N) * fctn(x_n);

end

approx = (b-a)/(3*N) * sum;    % compute the approximation from sum

end

function coeff = coefficient(n_input, N_input)
%  Compute the coefficient for the specified term according to
%  Simpson's Method.  In other words, return the correct coefficient
%  from this series of coefficients as used in Simpson's Method.
%  1,4,2,4,2,4, ... ,2,4,1
if n_input == 0 || n_input == N_input
    coeff = 1;
elseif rem(n_input, 2) == 0   % If n is even then coefficient = 2
    coeff = 2;
else
    coeff = 4;               % If n is odd then coefficient = 4
end

end

