微分变换法（DTM）求解FGM铁木辛柯梁固有频率时L/h=20不收敛的问题求助

最近几周我一直在尝试用**微分变换法（DTM）**求解功能梯度材料（FGM）铁木辛柯梁的固有频率。我开发了一段Python代码，通过牛顿-拉夫逊法求解特征方程来获取固有频率。目前当长厚比L/h=5时，结果收敛情况良好，但将长厚比调整为L/h=20时，就无法实现收敛了，想请各位帮忙分析问题所在，给点改进建议。

我的实现代码如下：

import sympy as sp
import numpy as np
from scipy import io, integrate, linalg, signal
from scipy.sparse.linalg import cg, eigs
import math
from math import *
from sympy import symbols, integrate, diff

# Define symbolic variables
d1, d2 = symbols('d1 d2 ')
z, x = sp.symbols('z x')
p, Ec, Em, nu = 0, 380., 70., 0.3
rhoc, rhom, alpha = 3960., 2702., 0.
l, Lh, s = 1., 5., 1.
h = l / Lh
b = l / s
ks = 5. / 6.

V1 = ((1./2.) + (z/h))**p;
Ez1 = (Ec - Em) * V1 + Em - (1/2) * alpha * (Ec + Em)
cnp = integrate(Ez1*z, (z, -h/2, h/2)) / (integrate(Ez1, (z, -h/2, h/2)));
print('c=', cnp)

V2 = ((1. / 2.) + ((z + cnp) / h))**p
# Define functions
Ez = (Ec - Em) * V2 + Em - (1/2) * alpha * (Ec + Em)
rhoz = (rhoc - rhom) * V2 + rhom - (1/2) * alpha * (rhoc + rhom)
fz = z
fz1 = diff(fz, z)
Q11 = Ez
Q55 = Ez / (2 * (1 + nu))
h1 = -h/2. - cnp; h2 = h/2. - cnp

# Perform symbolic integration
D11 = integrate(Q11 * z**2, (z, h1, h2))
A55a = ks * integrate(Q55 * fz1**2, (z, h1, h2))

I0 = integrate(rhoz, (z, h1, h2))
I2 = integrate(rhoz * z**2, (z, h1, h2))

# Compute coefficients
a11 = (-I0 * x) / A55a
a12 = -A55a / A55a
a21 = (-I2 * x + A55a) / D11
a22 = A55a / D11

print("Encastre-Encastre_(CC)")

# Define U and W as symbolic lists
Phi = [0 for _ in range(100)]  # Extend this as needed
W = [0 for _ in range(100)]  # Extend this as needed

# Initialize the given boundary conditions
Phi[0] = 0
Phi[1] = d1
W[0] = 0
W[1] = d2

for N in range(4, 26, 1):
    for K in range(0, N):
        W[K+2] = (a11 * W[K] + a12 * (K + 1) * Phi[K + 1]) / ((K+2)*(K+1))
        Phi[K+2] = (a21 * Phi[K] + a22 * (K + 1) * W[K + 1]) / ((K+2)*(K+1))
    
    # Summation equations
    f2 = sum(Phi[k] for k in range(N+1))
    f3 = sum(W[k] for k in range(N+1))
    
    f22 = sp.expand(f2)
    f33 = sp.expand(f3)

    # Collect terms involving d1, d2
    eqt2 = sp.collect(f22, {d1, d2})
    eqt3 = sp.collect(f33, {d1, d2})

    system = [eqt2, eqt3]

    # Extract the coefficients of d1, d2
    A = sp.Matrix([[eqt.coeff(d1) for eqt in system],
                   [eqt.coeff(d2) for eqt in system]])

    # Create the vector for the constants
    b = sp.Matrix([eqt.subs({d1: 0, d2: 0}) for eqt in system])

    # Compute the determinant of the matrix A
    det = A.det()
    ddet = diff(det, x)
    # Convert symbolic expressions to numerical functions
    f_numer = sp.lambdify(x, det, 'numpy')
    g_numer = sp.lambdify(x, ddet, 'numpy')

    # Define functions that return numerical values
    def f(x):
        return f_numer(x)  # Evaluate determinant at x

    def g(x):
        return g_numer(x)  # Evaluate derivative at x

    # Implementing Newton-Raphson Method
    def newtonRaphson(x0, ee, Nst):
        print('\n\n*** NEWTON RAPHSON METHOD ***')
        step = 1
        flag = 1
        condition = True
    
        while condition:
            if g(x0) == 0.0:  # Avoid division by zero
                print('Divide by zero error!')
                break
        
            x1 = x0 - f(x0) / g(x0)  # Now properly evaluates numerical values
        
            x0 = x1
            step += 1
        
            if step > Nst:
                flag = 0
                break
        
            condition = abs(f(x1)) > ee  # Ensure numerical comparison
    
        if flag == 1:
            print('\nRequired root is: %0.8f' % x1)
            Omega_bar = sqrt(x1) * l**2 * sqrt(rhom/Em) / h
            print(f'N={N:4.2f}   P={p:4.2f}   Lh={Lh:4.2f}   Omega_bar={Omega_bar}')
        else:
            print('\nNot Convergent.')

    # Input Section
    x0 = 0.009671  # input('Enter Guess: ')
    ee = 0.000001  # input('Tolerable Error: ')
    Nst = 1000     # input('Maximum Step: ')

    # Converting x0 and e to float
    x0 = float(x0)
    ee = float(ee)

    # Converting N to integer
    Nst = int(Nst)

    # Starting Newton Raphson Method
    newtonRaphson(x0, ee, Nst)

麻烦各位帮我分析下可能的原因，比如DTM的展开项数N的取值、牛顿-拉夫逊法的初始猜测值调整、FGM梁的力学建模细节，或者数值计算中的稳定性问题？有没有可行的改进方向能让L/h=20时也能收敛？

备注：内容来源于stack exchange，提问作者YouTldi