function y = rect_integrate(f, a, b)
N = 1001;
x = linspace(a,b,N);
delta = (b-a)/(N-1);
y = 0;
for i = 1 : N
    y = y + f(x(i)) * delta ;
end
end