function roots = find_roots(f, df, intervals, method, tol, imax)
    % FIND_ROOTS Refine roots using specified method
    % Inputs:
    %   f         - Function handle
    %   df        - Derivative function handle (for Newton's method)
    %   intervals - Nx2 array of bracketing intervals
    %   method    - 'bisection' or 'newton'
    %   tol       - Tolerance for convergence
    %   imax      - Maximum iterations
    % Outputs:
    %   roots - Array of refined root positions
    
    roots = [];
    for i = 1:size(intervals, 1)
        a = intervals(i, 1);
        b = intervals(i, 2);
        
        if strcmp(method, 'bisection')
            root = bisect(f, a, b, imax, tol);
        elseif strcmp(method, 'newton')
            x0 = (a + b) / 2; % Initial guess for Newton's method
            root = newton(x0, f, df, imax, tol);
        else
            error('Invalid method. Use ''bisection'' or ''newton''.');
        end
        roots = [roots; root];
    end
end