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