function roots = bisectionAllRoots(f, a, b, tol)
    % bisectionAllRoots: Finds all roots in the interval [a, b] using the
    % bisection method.
    %
    % Inputs:
    %   f   - function handle
    %   a   - start of interval
    %   b   - end of interval
    %   tol - tolerance for convergence
    %
    % Output:
    %   roots - array of roots found in the interval

    % Define the number of subintervals for root finding
    num_intervals = 1000; % Increase if more precision is needed
    x = linspace(a, b, num_intervals + 1);
    roots = [];

    % Iterate over subintervals
    for i = 1:num_intervals
        x1 = x(i);
        x2 = x(i + 1);

        % Check if there is a sign change in the interval
        if feval(f, x1) * feval(f, x2) < 0
            % Use bisection method to find the root in this subinterval
            root = bisect(f, x1, x2, 1000, tol); % Assuming bisect function is available
            roots = [roots, root];
        end
    end

    % Remove duplicates (due to numerical errors)
    roots = unique(roots);
end