function roots = BracketingMethod(f, d, imax, tol, interval, num_subintervals)
    % Inputs:
    % f: function handle
    % d: derivative function handle
    % imax: maximum number of iterations
    % tol: tolerance for convergence
    % interval: [a, b] search interval
    % num_subintervals: number of subintervals for bracketing

    % Output:
    % roots: array of roots found within the interval

    % Divide the interval into subintervals
    a = interval(1);
    b = interval(2);
    x_subintervals = linspace(a, b, num_subintervals+1);
    
    roots = []; % Initialize an array to store the roots
    
    % Loop through each subinterval
    for i = 1:num_subintervals
        x_left = x_subintervals(i);
        x_right = x_subintervals(i+1);
        
        % Check for a sign change in f(x) over the interval
        if feval(f, x_left) * feval(f, x_right) < 0
            % If a sign change is detected, use Newton's method
            
            % Initial guess is the midpoint of the interval
            x_guess = (x_left + x_right) / 2;
            counter = 1;
            abs_diff_between_guess = inf;

            % Newton's method iteration
            while (abs_diff_between_guess > tol) && (counter <= imax)
                f_at_guess = feval(f, x_guess);
                df_dx_at_guess = feval(d, x_guess);
                
                % Avoid division by zero
                if df_dx_at_guess == 0
                    break;
                end
                
                % Update guess
                x_next = x_guess - f_at_guess / df_dx_at_guess;
                abs_diff_between_guess = abs(x_guess - x_next);
                x_guess = x_next;
                counter = counter + 1;
            end
            
            % Check if the root is unique and within tolerance
            if abs(feval(f, x_guess)) < tol
                % Avoid duplicates
                if isempty(roots) || all(abs(roots - x_guess) > tol)
                    roots = [roots, x_guess]; % Add the root to the list
                end
            end
        end
    end
end
