多層導波路の数値解析に関する1992年の文献を読みました。この分野ではどうやら有名。だけど難解。オープンアクセスでないので詳しい内容を説明することはできませんが、手順を簡単にまとめたいと思います。多分このくらいの説明なら、この文献を引用するときに使っていても問題ないはず、という内容にとどめています。

転送行列法を用いて方程式を立てます
複素関数論の偏角の原理を用いて解の個数nを算出します
1の方程式をn次多項式で近似します
n次多項式を解きます

SEさんとかに見せたら怒られるコードだと思います。
層数が変わったら、関数M(gamma)を書き換える必要があります(笑)

DOI : https://ieeexplore.ieee.org/document/166774

import numpy as np
import cmath
from scipy.integrate import quad
ns,n1,n2,n3,n4,nc = 1.5,1.66,1.6,1.53,1.66,1.0
d1,d2,d3,d4 = 0.5,0.5,0.5,0.5 #um
wl = 0.6328 # um
k0 = 2*cmath.pi/wl

def M(gamma): # 1992 古い文献
    def Mi(ni,di):
        kap_i = cmath.sqrt((k0 * ni)**2 - gamma**2)
        return np.array([
            [cmath.cos(kap_i * di), 1j / kap_i * cmath.sin(kap_i * di)],
            [1j*kap_i * cmath.sin(kap_i * di), cmath.cos(kap_i * di)]
        ])
    return Mi(n1,d1) @ Mi(n2,d2) @ Mi(n3,d3) @ Mi(n4,d4)


def F(gamma):# 1992 古い文献
    # 行列を取得
    M_TMM = M(gamma)
    m11, m12 = M_TMM[0, 0], M_TMM[0, 1]
    m21, m22 = M_TMM[1, 0], M_TMM[1, 1]

    # 端部媒質の伝搬定数（複素数になる可能性あり）
    gamma_substrate = cmath.sqrt(gamma**2 - (k0 * ns)**2)
    gamma_cladding  = cmath.sqrt(gamma**2 - (k0 * nc)**2)

    # 関数の構築
    term1 = 1j*gamma_substrate * m11
    term2 = 1j*gamma_cladding  * m22
    term3 = -m21
    term4 = gamma_substrate * gamma_cladding * m12

    return term1 + term2 + term3 + term4



# 導関数を積分によって求める
def dF_Cauchy(gamma):
    def g(t,gamma):
        R = 1e-3

        return F(gamma+R*cmath.exp(1j*2*cmath.pi*t))/(R*cmath.exp(1j*2*cmath.pi*t))
    
    a = 0
    b = 1
    n = 100
    h = (b-a)/n # 微小区間
    G = h/2*((g(a,gamma)+g(b,gamma))+2*sum(g(a+h*i,gamma)for i in range(1,n)))
    
    return G

z1 = (1.501-1j)*k0
z2 = (2-1j)*k0
z3 = (2+1j)*k0
z4 = (1.501+1j)*k0



# 経路をパラメータ t ∈ [0, 1] で表現
def path(z_start, z_end, t):
    return z_start + (z_end - z_start) * t

def APM(P):
    def integrand_real(z_start, z_end):
        def f(t):
            gamma = path(z_start, z_end, t)
            dgamma_dt = z_end - z_start
            return (gamma**P * dF_Cauchy(gamma) / F(gamma)) * dgamma_dt
        return quad(lambda t: f(t).real, 0, 1)[0] + 1j * quad(lambda t: f(t).imag, 0, 1)[0]

    # 各辺の積分
    I_c1 = integrand_real(z1, z2)
    I_c2 = integrand_real(z2, z3)
    I_c3 = integrand_real(z3, z4)
    I_c4 = integrand_real(z4, z1)

    I_P = I_c1 + I_c2 + I_c3 + I_c4
    return I_P / (2 * cmath.pi * 1j)


sigma0 = APM(0)
print("解の個数は",sigma0,"個です。続いて解（実効伝搬定数ガンマ）を求めます。")

# 解の個数（次数）nの決定
n = int(round(sigma0.real))
sigma_list = []
for i in range(n+1):
    sigma_list.append(APM(i))


# ニュートンの恒等式 (べき乗和から多項式の係数を求める)
def newton_identities(n,p): 
    # 引数: 次数(n), べき乗和のリスト(p)
    e = [0]*(n+1)  # 対称式
    e[0] = 1
    poly = [1]  # 最高次の係数は1 
    for k in range(1,n+1):
        for i in range(1,k+1):
            e[k] += (-1)**(i-1) * e[k-i] * p[i]  # ニュートンの恒等式
        e[k] /= k
        poly.append((-1)**k * e[k])  
        
    return poly

poly = newton_identities(n,sigma_list)
S_z = np.poly1d(poly) # 多項式
print((S_z.r)/k0)

これの実部がTE0~TE3の実効屈折率です。文献とも合っているようです。
複素関数論の偏角の原理は、複素関数論の書籍でそんなに詳しく扱っていませんでした。演習問題で登場していました（同じ書籍で勉強した人もいるはず）。適切な経路で積分したら整数になって、それが解の個数になるって素晴らしいです。
ただ積分経路の取り方がかなり恣意的。答えがわかっているからこう設定できたようなもの。これを既存手法の課題として扱って、課題を解決したみたいな論文があったら知りたいです。

また時間があるときに書き足すかも。複素関数論って面白いですよね。