function roots = find_all_roots(f, d, a, b, imax, tol)
    % Input:
    %   f     - function handle for f(x)
    %   df    - function handle for derivative of f(x), f'(x)
    %   a, b  - interval limits
    %   imax  - maximum number of iterations
    %   tol   - tolerance for convergence
    % Output:
    %   roots - array of roots found in the interval
    
    % Initialize
    roots = []; % Variable to store all roots
    intervals = linspace(a, b, 1000); % Divide interval into small sections
    
    root_counter = 1; % Counter for dynamic variable naming
    
    for i = 1:length(intervals) - 1
        x1 = intervals(i);
        x2 = intervals(i+1);
        
        % Check for a sign change (indicating a root exists)
        if feval(f, x1) * feval(f, x2) < 0
            % Use Bisection to find an approximate root
            root_bisection = bisect(f, x1, x2, imax, tol);
            
            % Refine using Newton's method
            [root_newton, ~] = newton(f, d, root_bisection, imax, tol);
            
            % Add to roots array if unique
            if isempty(roots) || all(abs(roots - root_newton) > tol)
                roots = [roots; root_newton];
                
                % Dynamically assign variable for each root
                assignin('base', sprintf('root%d', root_counter), root_newton);
                root_counter = root_counter + 1;
            end
        end
    end
end

function [root, counter] = bisect(f, a, b, imax, tol)
    % Bisection Method Implementation
    counter = 0;
    x1 = a;
    x3 = b;
    f1 = feval(f, x1);
    f3 = feval(f, x3);
    
    if f1 * f3 > 0
        error('No root in the specified interval');
    end
    
    for counter = 1:imax
        x2 = (x1 + x3) / 2;
        f2 = feval(f, x2);
        if abs(f2) < tol || (x3 - x1) / 2 < tol
            root = x2;
            return;
        end
        if f1 * f2 < 0
            x3 = x2;
        else
            x1 = x2;
        end
    end
    root = (x1 + x3) / 2;
    warning('Maximum iterations reached, result may not be within tolerance');
end

function [root, counter] = newton(f, df, x_guess, imax, tol)
    % Newton's Method Implementation
    counter = 0;
    while counter < imax
        f_val = feval(f, x_guess);
        df_val = feval(df, x_guess);
        if abs(df_val) < tol
            error('Derivative near zero, method fails');
        end
        x_next = x_guess - f_val / df_val;
        if abs(x_next - x_guess) < tol
            root = x_next;
            return;
        end
        x_guess = x_next;
        counter = counter + 1;
    end
    root = x_guess;
    warning('Maximum iterations reached, result may not be within tolerance');
end