function I = simpson(f, a, b, n)
h = (b - a)/(2*n);
x = a:h:b;
y = f(x);
I = y(1) + y(2*n+1) + 4*(sum(y(2:2:(2*n)))) + 2*(sum(y(3:2:(2*n-1))));
I = I * h / 3;
end
