Python Error: Could not find root within given tolerance?

I'm trying to find the roots of an equation f that is the result of the determinant of a matrix. Arriving at this equation was simple, I generated the following code:

from sympy import *
from math import sqrt

def Rec_Coef(n, l, Q, rho, i):
    r1 = (1/2)*(1+sqrt(1-4*Q**2))
    r2 = (1/2)*(1-sqrt(1-4*Q**2))
    
    b = (rho*r1**2)/(r1-r2)
    
    if i == 1:
        q = (1/2)*(3+sqrt(9+16*Q**2*(l-1)*(l+2)))
    if i == 2: 
        q = (1/2)*(3-sqrt(9+16*Q**2*(l-1)*(l+2)))
        
    A = l*(l+1)    
        
    alpha = (n**2 + 2*(b+1)*n + 2*b + 1)*r1
    beta = (-2 + r2)*n**2 + (-2 - 2*b*(2-r2) - 4*rho*r1**2 + 6*r2)*n + (q - 2*(rho**2)*(r1**3) - 4*rho*r1**2*(1+b) - 2*b**2*(2-r1**(-1)) - r1*A - (3-2*b)*r2)
    gamma = (1+r2)*n**2 + (2*rho*r1*(1+2*r2)+2*b*(1+r2)-10*r2)*n + (rho*r1*(2 -12*r2 + (rho+2*b)*(1+2*r2)) -1 -q -2*b*(1+3*r2) + b**2*(16 + 8*r2 - (15 - 38*r2 + 26*r2**2)*r1**(-3)) + (A+13)*r2)
    delta = (-n**2 + 2*(3 - rho -b)*n - (9 + 4*rho*b - 6*b - 6*rho))*r2
    
    return alpha, beta, gamma, delta

N=4

α = []
β = []
γ = []
δ = []

l = 2 
Q = 0.2
i = 2

rho = symbols("𝜌")

for j in range(N):
    alpha, beta, gamma, delta = simplify(Rec_Coef(j, 2, 0.2, rho, 2))
    α.append(alpha)
    β.append(beta)
    γ.append(gamma)
    δ.append(delta)  

array = [[0]*(N+1) for row in range(N)]

array[0][0] = β[0]
array[0][1] = α[0]
array[1][0] = γ[1]
array[1][1] = β[1]
array[1][2] = α[1]

for i in range(2,N):
    array[i][i-1] = γ[i] - (array[i-1][i-1]/array[i-1][i-2])*δ[i]
    array[i][i] = β[i] - (array[i-1][i]/array[i-1][i-2])*δ[i]
    array[i][i+1] = α[i] 
    
matrix_aux = simplify(Matrix(array))

M_aux = [[0]*(N) for row in range(N)] 

for i in range(N):
    for j in range(N):
        M_aux[i][j]=(array[i][j])  

matrix = Matrix(M_aux)

f = matrix.det()

So I tried to find the roots in a few ways:

First

I used the mpmath function, turned my symbolic equation into a numerical equation, and then I guessed at the value. The reference value is (-0.17792, -0.74734*I).

import mpmath as mp

f = lambdify(rho, matrix.det())

mp.findroot(f, -0.17-0.75j)

However, I got the following error:

ValueError: Could not find root within given tolerance. (2.01839854512653723801e-14 > 2.16840434497100886801e-19) Try another starting point or tweak arguments.

Second

I used the solve function

f = matrix.det()
sol = solve(f, rho)
sol = 
展开收缩
for s in sol:
    print(s)

But this time, even with a low number of iterations used at the beginning of the code (N=4), the cell keeps loading without stopping and does not generate any result or code.

How can I solve this problem and generate the complex roots easily?

Let me start by saying that I have no idea what problem you are trying to solve, hence I'm going to assume that your equation f is built correctly.
Also, I'm unable to solve your problem, but hopefully you'll gain some insight about why your attempts failed.

To gain more intuition about the function, I plotted its absolute value in the complex plane using the sympy plot backends module:

import numpy as np
plot_complex(f, (rho, -2-2j, 2+2j), threed=True, backend=KB,
             tz=lambda t: np.log(t)/10, n=500)

Note that the above command evaluates the function over 500*500=250000 points. Watch out for memory consumption! In the above picture, every spike towards the bottom is a zero.

Looking closer in the following region, and increasing a lot the resolution of the plot (3000*3000=9M evaluation points, don't run this command as it consumes over 15GB of memory):

plot_complex(f, (rho, -0.1-1j, -0.35-0.75j), threed=True, backend=KB,
             tz=lambda t: np.log(t)/30, n=3000)

See that noise? Let's zoom in:

In this region, it appears that your functions has "infinitely many zeros" and "infinitely many saddle points". This is a nightmare situation: no matter your initial condition, I believe it's impossible to find a zero with any numerical algorithm, let alone symbolic algorithms like solve or solveset. The distance between zeros is sooo small that I don't see how an algorithm can properly resolve the situation... Maybe findroot has some option to specify the tolerance???

By the way, sympy's nsolve is a wrapper to mpmath's findroot, so you can avoid the lambdification step!

Now, let's zoom into a different region:

plot_complex(log(f), (rho, -0.65-0.95j, -0.6-0.9j), threed=True, backend=KB,
             tz=lambda t: np.log(t)/100, n=3000)

See that noise? Let's zoom in a little bit:

It appears that this is a region with infinitely many zeros, hence infinitely many saddle points... Again, a nightmare...