How to find complex roots of an equation without analytical solution?

I have a complex function f which is given in terms of functions:

from sympy import *

def alpha(n, rho, l):
    return n**2 + (2*rho + 2)*n + 2*rho + 1

def beta(n, rho, l):
    return -(2*n**2 + (8*rho + 2)*n + 8*rho**2 + 4*rho + l*(l+1) -3)

def gamma(n, rho, l):
    return n**2 + 4*rho*n + 4*rho**2 - 4

With these functions I discretize a certain continued fraction:

def ContinuedFrac(n, rho, l):
    cf = 0
    for i in reversed(range(1,n+1)):
        cf = cf + beta(i, rho, l)/alpha(i, rho, l)
        cf = (-gamma(i, rho, l)/alpha(i, rho, l)) / cf
    return cf

And finally, I assemble the equation to be solved:

omega = symbols("ω")
n = 5
rho_guess = -I*omega
l=2

f = ContinuedFrac(n, rho_guess, l) + beta(0, rho_guess, l)/alpha(0, I*omega, l)

My goal is to solve the function f for the variable omega in order to find the complex roots, for example (0.7473, -0.1779). The problem is that this function does not apparently have an analytic solution, but has infinite numerical roots.

I tried to solve it using the lambdfy function, but with them, I need to give an initial value and I would like to know if it is possible to directly find a list with several omega roots of the function f.

Below is what I tried to do. It works, but I only find a function according to a certain initial value.

import mpmath as mp

omega = symbols("ω")

def findomega(n, rho_guess, l):
    return lambdify(omega, ContinuedFrac(n, rho_guess, l) + beta(0, rho_guess, l)/alpha(0, rho_guess , l))

f = findomega(15, -I*omega, 3)
mp.findroot(f, 1.16-0.56j)

output: mpc(real='1.1652715765041771', imag='-0.56260091739349236')

> ... and I would like to know if it is possible to directly find a list with several omega roots of the function f.

Not directly, but with a bit of effort you could get close to it. Here is how I would do it...

I'm going to use the Sympy Plotting Backend module, in particular the plot_complex function to create some kind of domain coloring plot. For this particular case, I'm choosing a black and white color scheme highlighting the zeros (black) and poles (white).

from spb import plot_complex

n = 5
l = 2
f = ContinuedFrac(n, rho_guess, l) + beta(0, rho_guess, l)/alpha(0, I*omega, l)
plot_complex(
    f, (omega, -2-2j, 2+2j),
    grid=False, coloring="k", n=500
)

So, black spots are zeros, white spots are poles (in general, white spots are zones of higher magnitude. In this particular case, it happens that they are poles).
From this picture you can easily extract initial guesses for the zeros, and use sympy's nsolve to find the numerical root. For example:

nsolve(f, omega, 0.54-0.1j)
# out: 0.533191990477995−0.100182203167779𝑖

Using this particular color scheme, I believe one could easily write a function to locate the darkest spots on the image, in order to automatically extract initial guesses and provide them to nsolve. For example:

p = plot_complex(...
xx, yy, _, _, img, _ = p[0].get_data()

# your procedure to extract zeros from img. For example:
# 1. `img` is a black and white RGB array (n x m x 3). Convert it to
#    a grey scale array (n x m), so each element gives the grey value.
#    0=black, 1=white, ...
# 2. Increase the contrast of the image (you can use PIL): this
#    makes the light-grey areas become white, and dark grey areas
#    become black.
# 3. Run a clustering algorithm on black levels: this will find
#    the centers of the black spots in image coordinates.
# 4. Convert the centers in data coordinates (you can use the x
#    and y extent of the image).
# 5. Use the centers in data coordinates as initial guesses
#    for nsolve.

That's for your to try

Since it appears that you are interested in studying the function by varying n and l, I would also use the interactive capabilities of the plotting module to create a sliding window so that I can easily move/zoom around the interested region of the plot.

If you are somehow able to rewrite your function using SymPy's Sum, it would also be possible to use symbolic n and l, declare them as parameters in the plotting function, in order to quickly change their values with a widget (for example, ipywidgets.BoundedIntText which would render as a spinner for integer values).

%matplotlib widget
from spb import plot_complex, prange

n = 5
l = 2
f = ContinuedFrac(n, rho_guess, l) + beta(0, rho_guess, l)/alpha(0, I*omega, l)

# sliding window of the plot: center x, center y, width, height
cx, cy, w, h = symbols("cx, cy, w, h")

plot_complex(
    f, prange(omega, (cx - w / 2) + (cy - h / 2)*I, (cx + w / 2) + (cy + h / 2)*I),
    grid=False, coloring="k", n=500,
    params={
        cx: (0, -5, 5),
        cy: (0, -5, 5),
        w: (4, 0, 10),
        h: (4, 0, 10),
    }
)