I first did a numerical calculation with numpy and the result was correct, but due to insufficient precision, it could lead to errors. So I used mpmath for the calculation instead, but I found that the result was different.

import sympy as sy
import numpy as np
from sympy import Poly
from functools import lru_cache

y, omega = sy.symbols('y omega')

def lambda0_f(y):
    return 2 * y * (1 - 1j * omega) / (1 - y**2)

def s0_f(y):
    return (1 - 2 * 1j * omega - 2 * omega**2) / (2 * (1 - y**2))

y0_val = 0.4
num_val = 3

lambda0_series = sy.series(lambda0_f(y), y, y0_val, num_val+1).removeO()
s0_series = sy.series(s0_f(y), y, y0_val, num_val+1).removeO()

def c0i(i0):
    if i0 == 0:
        return lambda0_f(y0_val)
    else:
        return lambda0_series.coeff((y-y0_val)**i0) 

def d0i(i0):
    if i0 == 0:
        return s0_f(y0_val)
    else:
        return s0_series.coeff((y-y0_val)**i0)
  
c0i_coeffs = [Poly(c0i(i), omega).all_coeffs() for i in range(num_val+1)]
d0i_coeffs = [Poly(d0i(i), omega).all_coeffs() for i in range(num_val+1)]

c0i_coeffs_np = [np.array(coef, dtype=np.complex128) for coef in c0i_coeffs]
d0i_coeffs_np = [np.array(coef, dtype=np.complex128) for coef in d0i_coeffs]

def convolve_polys(poly1, poly2):
    return np.convolve(poly1, poly2)

def poly_add(poly1, poly2):
    len_diff = len(poly1) - len(poly2)
    if len_diff > 0:
        poly2 = np.pad(poly2, (len_diff, 0), mode="constant")
    elif len_diff < 0:
        poly1 = np.pad(poly1, (-len_diff, 0), mode="constant")
    return poly1 + poly2

@lru_cache(maxsize=None)
def compute_cni_dni_coeffs(i0, num, c0_coeffs, d0_coeffs):
    if num == 0:
        return c0_coeffs[i0], d0_coeffs[i0]

    c_result = np.zeros(num, dtype=np.complex128)
    d_result = np.zeros(num, dtype=np.complex128)

    for k in range(i0 + 1):
        c_result, d_result=poly_add(c_result, convolve_polys(np.array(c0_coeffs[k]), np.array(compute_cni_dni_coeffs(i0 - k, num - 1, tuple(map(tuple, c0_coeffs)), tuple(map(tuple, d0_coeffs)))[0]))),poly_add(d_result, convolve_polys(np.array(d0_coeffs[k]), np.array(compute_cni_dni_coeffs(i0 - k, num - 1, tuple(map(tuple, c0_coeffs)), tuple(map(tuple, d0_coeffs)))[0])))

    if i0 + 1 <= num_val:
        c_result, d_result=poly_add(c_result, poly_add((i0 + 1) * np.array(compute_cni_dni_coeffs(i0 + 1, num - 1, tuple(map(tuple, c0_coeffs)), tuple(map(tuple, d0_coeffs)))[0]), np.array(compute_cni_dni_coeffs(i0, num - 1, tuple(map(tuple, c0_coeffs)), tuple(map(tuple, d0_coeffs)))[1]))),poly_add(d_result, (i0 + 1) * np.array(compute_cni_dni_coeffs(i0 + 1, num - 1, tuple(map(tuple, c0_coeffs)), tuple(map(tuple, d0_coeffs)))[1]))


    return tuple(c_result), tuple(d_result)

cni_coeffs, dni_coeffs = compute_cni_dni_coeffs(0, num_val, tuple(map(tuple, c0i_coeffs_np)), tuple(map(tuple, d0i_coeffs_np)))

cni_coeffs, dni_coeffs = compute_cni_dni_coeffs(0, num_val, tuple(map(tuple, c0i_coeffs_np)), tuple(map(tuple, d0i_coeffs_np)))

print("cni_coeffs:", cni_coeffs)
print("dni_coeffs:", dni_coeffs)

from functools import lru_cache
import mpmath
import numpy as np
from mpmath import mp, mpc
import sympy as sy
from sympy import symbols, Poly

ndps=20
mp.dps = ndps 

y, rh, r2, omega, l ,M,Q= sy.symbols('y r_h r_2 omega ell M Q')

def lambda0_f(y):
    return 2 * y * (1 - 1j * omega) / (1 - y**2 + 1e-15)


def s0_f(y):
    return (1 - 2 * 1j * omega - 2 * omega**2) / (2 * (1 - y**2 + 1e-15))

y0_val = 0.4
num_val =3

lambda0_series = sy.series(lambda0_f(y), y, y0_val, num_val+1).removeO()
s0_series = sy.series(s0_f(y), y, y0_val, num_val+1).removeO()

def c0i(i0):
    if i0 == 0:
        return lambda0_f(y0_val)
    else:
        return lambda0_series.coeff((y-y0_val)**i0) 


def d0i(i0):
    if i0 == 0:
        return s0_f(y0_val)
    else:
        return s0_series.coeff((y-y0_val)**i0)

c0i_coeffs = [Poly(c0i(i),omega).all_coeffs() for i in range(num_val+1)]
d0i_coeffs = [Poly(d0i(i),omega).all_coeffs() for i in range(num_val+1)]

c0i_coeffs_np = [np.array(coef, dtype=np.complex128) for coef in c0i_coeffs]
d0i_coeffs_np = [np.array(coef, dtype=np.complex128) for coef in d0i_coeffs]


def numpy_to_mpmath(x):
    return [[mpc(np.real(y),np.imag(y)) for y in z ] for z in x ]

lambda_coes_mpmath = [[mpc(np.real(y),np.imag(y)) for y in z ] for z in c0i_coeffs_np ]
s_coes_mpmath = [[mpc(np.real(y),np.imag(y)) for y in z ] for z in d0i_coeffs_np ]

def poly_add(poly1, poly2):
    if len(poly2) > len(poly1):
        poly1, poly2 = poly2, poly1
    result = [mpc(0) for _ in range(len(poly1))]
    for i in range(len(poly2)):
        result[i] = poly1[i] + poly2[i]
    for i in range(len(poly2), len(poly1)):
        result[i] = poly1[i]
    return result

def convolve_polys(p1, p2):
    result_length = len(p1) + len(p2) - 1
    result = [mpc(0) for _ in range(result_length)]
    for i in range(len(p1)):
        for j in range(len(p2)):
            result[i + j] += p1[i] * p2[j]
    return result

@lru_cache(maxsize=None)
def compute_cni_dni_coeffs(i0, num, c0_coeffs, d0_coeffs):
    if num == 0:
        return c0_coeffs[i0], d0_coeffs[i0]

    c_result = [mpc(0) for _ in range(num)]
    d_result = [mpc(0) for _ in range(num)]

    for k in range(i0 + 1):
        c_k, d_k = compute_cni_dni_coeffs(i0 - k, num - 1, c0_coeffs, d0_coeffs)
        c_result,d_result = poly_add(c_result, convolve_polys(c0_coeffs[k], c_k)),poly_add(d_result, convolve_polys(d0_coeffs[k], c_k))
    if i0 + 1 <= num_val:
        c_next, d_next = compute_cni_dni_coeffs(i0 + 1, num - 1, c0_coeffs, d0_coeffs)
        _,d_c0=compute_cni_dni_coeffs(i0, num - 1, c0_coeffs, d0_coeffs)
        c_result = poly_add(c_result, poly_add(convolve_polys(c_next, [mpc(i0 + 1,0)]), d_c0))
        d_result = poly_add(d_result, convolve_polys(d_next, [mpc(i0 + 1,0)]))

    return c_result, d_result

cni_coeffs, dni_coeffs = compute_cni_dni_coeffs(0, num_val, tuple(map(tuple, lambda_coes_mpmath)), tuple(map(tuple, s_coes_mpmath)))

Please help me get the right results from mpmath, thanks!