% Main script to evaluate the integral for f(x) = 1/x
a = 1; % lower limit
b = 2; % upper limit
intervals = [5, 50, 100]; % different numbers of intervals
results = zeros(length(intervals), 2); % To store results and errors

for i = 1:length(intervals)
    n = intervals(i);
    I_numeric = composite_trap('f_reciprocal', a, b, n);
    results(i, 1) = I_numeric; % Store numerical result

    % Calculate the analytic value
    analytic_value = log(2); % ln(2)
    results(i, 2) = analytic_value; % Store analytic result

    % Calculate percentage error
    percentage_error = abs((I_numeric - analytic_value) / analytic_value) * 100;
    disp(['Numerical integral for ', num2str(n), ' intervals: ', num2str(I_numeric)]);
    disp(['Analytic integral: ', num2str(analytic_value)]);
    disp(['Percentage error: ', num2str(percentage_error), '%']);
end
% Create a table for results
T = table(intervals', results(:, 1), results(:, 2), ...
          'VariableNames', {'Intervals', 'Numerical Integral', 'Analytic Integral'});
disp(T);

% Determine the number of sub-intervals needed for 6 decimal places of accuracy
tolerance = 1e-6; 
n = 1; 
error = Inf; 

while error > tolerance
    I_numeric = composite_trap('f_reciprocal', a, b, n);
    analytic_value = log(2);
    error = abs((I_numeric - analytic_value) / analytic_value); 
    n = n + 1; 
end

disp(['Number of sub-intervals needed for 6 decimal places of accuracy: ', num2str(n-1)]);