function y = simp_integrate(f, a, b)
N = 1001; % N is odd for Simpson's method if starting index is 1
x = linspace(a,b,N);
delta = (b-a)/(N-1);
y = 0;
y = y + f(x(1)) + f(x(N));
for i = 2 : N-1
    if (mod(i,2) == 0)
    y = y + 4*f(x(i)) ;
    else
    y = y + 2*f(x(i)) ;
    end
end
y = y * delta / 3;
end