"Bisection method: Finds a zero in the interval [a, b]"
function bisection(f::Function, a::Float64, b::Float64; tol::Float64 = 1e-8, max_iter::Int = 1000)::Float64
    fa::Float64 = f(a)
    fb::Float64 = f(b)
    
    # Ensure the endpoints bracket a root.
    if fa * fb > 0.0
        error("No sign change detected in the interval: f(a)*f(b) > 0")
    end

    c::Float64 = a
    for _ in 1:max_iter
        c = (a + b) / 2.0
        fc::Float64 = f(c)
        
        # Check for convergence using the function value or the interval width.
        if abs(fc) < tol || abs(b - a) < tol
            return c
        elseif fa * fc < 0.0
            b = c
            fb = fc
        else
            a = c
            fa = fc
        end
    end
    return c  # Return the current approximation if max iterations reached.
end

"Function to find all zeros in [x_start, x_end] by scanning for sign changes."
function find_all_zeros(f::Function, x_start::Float64, x_end::Float64; partitions::Int = 1000, tol::Float64 = 1e-8)::Vector{Float64}
    zeros_list::Vector{Float64} = Float64[]
    Δ::Float64 = (x_end - x_start) / partitions
    x_prev::Float64 = x_start
    f_prev::Float64 = f(x_prev)
    
    for i in 1:partitions
        x_curr::Float64 = x_start + i * Δ
        f_curr::Float64 = f(x_curr)
        
        # Check for a sign change or an exact zero.
        if f_prev * f_curr < 0.0 || f_prev == 0.0
            try
                root::Float64 = bisection(f, x_prev, x_curr; tol=tol)
                # Add the root if it's new (avoid duplicates from neighboring intervals)
                if isempty(zeros_list) || abs(root - zeros_list[end]) > tol
                    push!(zeros_list, root)
                end
            catch err
                # Skip this interval if bisection fails (should not occur if a sign change exists)
            end
        elseif f_curr == 0.0
            if isempty(zeros_list) || abs(x_curr - zeros_list[end]) > tol
                push!(zeros_list, x_curr)
            end
        end
        
        x_prev = x_curr
        f_prev = f_curr
    end
    return zeros_list
end