function [I] = mc2(f,a,b,N)
X = a + (b-a).*rand(N,1);
Y = a + (b-a).*rand(N,1);
under_graph = find(Y<=f(X));
I = (length(under_graph)/N) * (b-a) * (b - a);


end