% Main script to evaluate the integral for f(x) = 1/x
a = 1; 
b = 2; 
intervals = [5, 50, 100]; 

results = zeros(length(intervals), 2); 

for i = 1:length(intervals)
    n = intervals(i);
    % Ensure n is even for Simpson's Rule
    if mod(n, 2) ~= 0
        n = n + 1;
    end
    I_numeric = simpson_rule('f_reciprocal', a, b, n);
    results(i, 1) = I_numeric;

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

    % 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 = 2; 
error = Inf; 

% Initialize the analytic value once before the loop
analytic_value = log(2);

% Loop until the error is within the specified tolerance
while error > tolerance
    I_numeric = simpson_rule('f_reciprocal', a, b, n); 
    error = abs((I_numeric - analytic_value) / analytic_value); 
    n = n + 2; 
end

% Display the number of intervals needed to achieve the desired accuracy
disp(['Number of sub-intervals needed for 6 decimal places of accuracy: ', num2str(n-2)]);