function [ I ] = MySimpson( f, a, b, n )
% Simspon formula evaluated
% n number of sections of length 2h with h = x_i+1 - x_i

if a > b 
    temp = a;
    a = b;
    b = temp;
    disp('a and b interchanged')
end
    
h = (b-a)/(2*n)

% Special cases
I = 1/3*(f(a) + f(b))*h + 4/3 * f(b -h);

for i = 1:n-1
    % inner points (all but last one)
    I = I +  4/3 * f(a + ( 2*i-1) * h);
    % touching points 
    I = I +  2/3 * f(a +   2*i    * h);
end

I = I * h

end