Решение нелинейной системы ОДУ с краевыми условиями на Python

Я пытаюсь решить эту систему дифференциальных уравнений:

Функции должны вести себя следующим образом, когда r-> бесконечность (F = g = lambda = const.):

Недавно у меня была похожая проблема, и один пользователь любезно предоставил мне алгоритм, с помощью которого я мог решить это ОДС. Теперь я попытался применить алгоритм TDM к новой задаче, но получил только тривиальное решение. Кто-нибудь знает, что я сделал не так, или как я могу решить эту проблему?

import numpy as np
import matplotlib.pyplot as plt

# Функция TDMA для решения уравнения трёхдиагональной матрицы
def tdma(a, b, c, d):
    n = len(d)
    p = np.zeros(n)
    q = np.zeros(n)
    x = np.zeros(n)

    # Прямой проход
    i = 0
    denominator = b[i]
    p[i] = -c[i] / denominator
    q[i] = d[i] / denominator
    for i in range(1, n):
        denominator = b[i] + a[i] * p[i - 1]
        if abs(denominator) < 1.0e-10:
            print("Нет решения")
            return x
        p[i] = -c[i] / denominator
        q[i] = (d[i] - a[i] * q[i - 1]) / denominator

    # Обратный проход
    x[n - 1] = q[n - 1]
    for i in range(n - 2, -1, -1):
        x[i] = p[i] * x[i + 1] + q[i]

    return x

# Параметры
N = 1000  # Количество ячеек (1-N); добавлены границы 0 и N+1
rmin, rmax = 0.0, 100.0  # Минимальное и максимальное r
dr = rmax / N  # Ширина ячейки
lamda = 10
g = 1.0  # Вы можете настроить этот параметр, если необходимо

# Граничные условия на F и W
FL, FR = 0, 1
WL, WR = 0, 1 / rmax

# Расположение ячеек / узлов
r = np.linspace(rmin - 0.5 * dr, rmax + 0.5 * dr, N + 2)
r[0], r[N + 1] = rmin, rmax
rv = np.linspace(rmin, rmax, N + 1)

# Коэффициенты
awF = np.zeros(N + 2)
aeF = np.zeros(N + 2)
awW = np.zeros(N + 2)
aeW = np.zeros(N + 2)

# Основные ячейки
aeF[1:N + 1] = 1 / dr ** 2
aeW[1:N + 1] = rv[1:N + 1] ** 2 / dr ** 2
awF[1:N + 1] = aeF[0:N]
awW[1:N + 1] = aeW[0:N]

# Границы
awF[1], awW[1], aeF[N], aeW[N] = 0, 0, 0, 0
awFL = 2 / dr ** 2
aeFR = 2 / dr ** 2
awWL = 2 * rmin ** 2 / dr ** 2
aeWR = 2 * rmax ** 2 / dr ** 2

# Начальные значения
F = np.ones(N + 2) * 1e-5
W = np.ones(N + 2) * 1e-5

# Граничные условия
F[0], F[N + 1] = FL, FR
W[0], W[N + 1] = WL, WR

# Фактор недоразмещения
alpha = 0.5
niter = 0
s = np.zeros(N + 2)
sp = np.zeros(N + 2)

# Цикл итераций
for _ in range(2000):
    niter += 1
    F_target = F.copy()
    W_target = W.copy()

    ######################
    # Обновление уравнения W #
    ######################
    for i in range(1, N + 1):
        K = 1 - g * r[i] * W[i]
        H = g * r[i] * F[i]

        # Ограничение значений для предотвращения переполнения
        K = np.clip(K, -1e5, 1e5)
        H = np.clip(H, -1e5, 1e5)

        s[i] = K * (K**2 - 1) + H**2 * K
        sp[i] = -2 * (1 - g * r[i]) / r[i]**2

        if np.isnan(K) or np.isnan(H):
            print(f"Числовая проблема при r = {r[i]}, K = {K}, H = {H}")
            break

    s[1] += awFL * FL
    sp[1] -= awFL
    s[N] += aeFR * FR
    sp[N] -= aeFR

    ap = awF + aeF - sp
    F_target[1:N + 1] = tdma(-awF[1:N + 1], ap[1:N + 1], -aeF[1:N + 1], s[1:N + 1])

    ######################
    # Обновление уравнения F #
    ######################
    for i in range(1, N + 1):
        K = 1 - g * r[i] * W[i]
        H = g * r[i] * F[i]

        # Ограничение значений для предотвращения переполнения
        K = np.clip(K, -1e5, 1e5)
        H = np.clip(H, -1e5, 1e5)

        s[i] = 2 * H * K**2 + lamda * H * ((H**2 / g**2) - r[i]**2 * F[i]**2)

        if np.isnan(K) or np.isnan(H):
            print(f"Числовая проблема при r = {r[i]}, K = {K}, H = {H}")
            break

    s[1] += awWL * WL
    sp[1] -= awWL
    s[N] += aeWR * WR
    sp[N] -= aeWR

    ap = awW + aeW - sp
    W_target[1:N + 1] = tdma(-awW[1:N + 1], ap[1:N + 1], -aeW[1:N + 1], s[1:N + 1])

    # Вычислить изменение
    change = (np.linalg.norm(F_target - F) + np.linalg.norm(W_target - W)) / N

    # Недоразмещение обновления
    F = alpha * F_target + (1 - alpha) * F
    W = alpha * W_target + (1 - alpha) * W

    if change < 1.0e-10:
        break

print(niter, " итераций ")

# Построить результаты
plt.plot(r, F, label="F(r)", color="g")
plt.plot(r, W, label="W(r)", color="r")
plt.legend(loc="upper right")
plt.xlabel("r")
plt.ylim(0)
plt.xlim(0, 20)
plt.grid(True)
plt.show()