function Integral

N = 10;
x_axis = zeros(N,1);
x_0 = 1; x_1 = 2; dx = (x_1-x_0)/N;
f = zeros(N+1,1);
for ix = 1:N+1
    x_j = x_0+dx*(ix-1);
    x_axis(ix,1) = x_j;
    f(ix,1) = x_j^(-1/2)*exp(-x_j);
end

% plot(x_axis,f)
% stop

int_f = 0;
for ix=1:N
    int_f = int_f + dx*(f(ix,1)+f(ix+1,1))/2;
end

error = 1;
while error > 10^(-5)
    N = N+1;

    clear x_axis f
    x_axis = zeros(N,1);
    x_0 = 1; x_1 = 2; dx = (x_1-x_0)/N;
    f = zeros(N+1,1);
    for ix = 1:N+1
        x_j = x_0+dx*(ix-1);
        x_axis(ix,1) = x_j;
        f(ix,1) = x_j^(-1/2)*exp(-x_j);
    end

    int_f_Np1 = 0;
    for ix=1:N
        int_f_Np1 = int_f_Np1 + dx*(f(ix,1)+f(ix+1,1))/2; %dx*f(ix,1); % dx*(f(ix,1)+4*f(ix+1,1)+f(ix+2,1))/3; %
    end
    
    error = abs(int_f_Np1-int_f);
    
    int_f = int_f_Np1;
    
end

N