function [x, dx, il, ir, I_num]=indef_integral(a, b, N, g)

x=linspace(a,b,N); dx=x(2)-x(1);
il=(1:1:N-1); ir=(2:1:N);
I_num(1)=0
trap=dx*(g(x(il))+g(x(ir)))/2;
I_num(ir)=cumsum(trap);

end