Мне нужно получить такое же решение, как представл...

Да. Я проверил основной код до строки ClearAll[scanResidual]; сам scanResidual и всё, что идёт после него, я не анализировал, как вы просили.

Главное: проблема не одна. В коде есть несколько независимых ошибок, и две из них сами по себе достаточны, чтобы shoot возвращал $Failed.

Что ломает код

1. Самая серьёзная программная ошибка — смешаны функции и их численные значения.

Например, у вас определено примерно так:

wl
nuPrimeFun[r_, phi_, phip_, lam_] := ...
    Exp[lam[r]]
    phip[r]^2
    ...

а вызывается:

wl
nuPrimeFun[r, phi[r], F[r], lam[r]]

После подстановки Mathematica фактически получает конструкции вида

wl
Exp[lam[r][r]]
F[r][r]^2
phi[r][r]

То есть аргумент lam_ уже равен значению lam[r], но затем вы второй раз применяете к нему [r].

Эта ошибка была сразу в нескольких местах:

wl
Aexpr
nuPrimeFun
lambdaPrimeFun
C33
rhsVals

Именно поэтому правая часть системы перестаёт быть численной для NDSolve.

В исправленной версии везде разделены численные аргументы:

wl
Aexpr[rr_, ph_, pp_, la_] := ...
Exp[la]
ph
pp

а не la[rr], ph[rr], pp[rr].

---

2. shoot[0.2] для ваших параметров вообще не должен работать.

В статье $\phi_S$ — это именно параметр стрельбы на горизонте, который подбирается так, чтобы на бесконечности получилось $\phi\to v$. Авторы прямо это пишут. (2207.03944v2.pdf)

То есть

wl
phiS = v

ставить нельзя.

Для ваших параметров

wl
ell = 3;
mu = 0.1;
h = 0.5;
v = 0.2;

режим как раз соответствует фиолетовой кривой Fig. 3, $\mu=0.1,\ell=3,v=0.2$. (2207.03944v2.pdf)

В этом режиме поле около горизонта почти равно нулю, а не $v=0.2$.

Я отдельно подставил ваши исходные значения rS=1.01, ell=3, phiS=0.2 в (30)+(32). Для phiS=0.2 имеем $V=V_{,\phi}=0$, и уравнение (32) даёт для $\phi'_S$

$$\phi'_S \simeq 0.03507\pm1.71454\,i .$$

Реального $\phi'_S$ вообще нет. Поэтому ваш

wl
FindRoot[..., {x, 0.1}]

и выдаёт FindRoot::lstol. Это не проблема начального приближения 0.1: при phiS=0.2 искомого вещественного корня нет.

По независимой численной проверке исправленной системы для параметров Fig. 3 и rmax=50 физический параметр стрельбы находится примерно около

$$\boxed{\phi_S \approx 6.4\times10^{-9}},$$

а не около 0.2 и тем более не около 0.87.

---

3. phiPrimeS нельзя надёжно определять одним FindRoot[{x,0.1}].

Уравнение (32), совместно с (30), является алгебраическим по $\phi'_S$. В зависимости от $\phi_S$ оно может иметь:

• два вещественных корня;
• один предельный корень;
• вообще не иметь вещественных корней.

Причём нужно ещё выбрать физическую ветвь с

$$b_1>0.$$

В статье именно (30) задаёт $b_1$, а (32) — граничное условие для производной поля. (2207.03944v2.pdf) (2207.03944v2.pdf)

Я заменил этот FindRoot на NSolve[..., Reals], отбрасываю корни с b1<=0 и выбираю физическую ветвь, непрерывно связанную с приближением

wl
24 ell^2 phiS/rS^3 + rS Vp[phiS]

---

4. rS = 1.01 неверно, если вы хотите воспроизвести графики статьи.

Авторы переходят к безразмерной координате

$$r\to r/r_S,$$

следовательно, после этого горизонт находится в

$$r_S=1.$$

(2207.03944v2.pdf)

У вас одновременно стояло

wl
rS = 1.01;
eps = 0.01;

а затем интегрирование начиналось с

wl
rS + eps

то есть фактически с 1.02.

Нужно:

wl
rS = 1;
eps = 1/100;
r0 = rS + eps;   (* 1.01 *)

eps и есть отступ от горизонта.

---

5. У вас закомментирована обязательная проверка b1 > 0.

Было:

wl
(* If[bb <= 0, Return[$Failed]]; *)

Но затем используется

wl
lam[r0] == -Log[bb eps]

Если bb<=0, решение немедленно становится комплексным. Проверку нужно вернуть.

---

6. Quiet@Check[..., $Failed] скрывает настоящую причину ошибки.

Конструкция

wl
Quiet @ Check[
    NDSolve[...],
    $Failed
]

очень неудобна для такого расчёта: Check может превратить вычисление в $Failed при сообщении от NDSolve, а Quiet одновременно скрывает это сообщение.

Поэтому вы видите только $Failed и не видите, что произошло внутри.

Я это убрал: теперь NDSolveValue показывает настоящие диагностические сообщения, а $Failed возвращается только если действительно не получилось получить четыре InterpolatingFunction.

---

7. MachinePrecision здесь недостаточно надёжен.

У вас:

wl
mu = 0.1;
h = 0.5;
rS = 1.01;

то есть машинные числа, а параметр стрельбы для рассматриваемого решения порядка $10^{-9}$. Кроме того, наружная стрельба является очень чувствительной separatrix-задачей: надо подавить растущую асимптотическую моду.

И сам Mathematica у вас сообщает:

text
You may need more than MachinePrecision digits ...

Я перевёл параметры в точные числа:

wl
mu = 1/10;
v  = 1/5;
h  = mu/v;

и использовал

wl
WorkingPrecision -> 50
AccuracyGoal -> 28
PrecisionGoal -> 22

---

8. В findPhiS у вас неправильное уравнение.

Вы определяете

wl
residual[phiS] = phi[rmax] - v

но затем ищете

wl
residual[phiS] == 0.1

Для условия (3)

$$\phi_\infty=v$$

должно быть строго

wl
residual[phiS] == 0

Условие $\phi_\infty=v$ дано в статье в (3), а метод стрельбы именно подбирает $\phi_S$, чтобы его выполнить. (2207.03944v2.pdf) (2207.03944v2.pdf)

Кроме того, стартовать FindRoot с

wl
phi0 = 0.0

для этой задачи неудобно. Для ваших параметров корень порядка $10^{-9}$, поэтому я использовал две стартовые точки:

wl
{x, 6*10^-9, 7*10^-9}

---

9. Ошибки есть и в ваших тестовых Plot до scanResidual.

Вот это неверно:

wl
Plot[nu /. shoot[0.87], ...]

Во-первых, надо вычислять nu[rr], а не символ nu. Во-вторых, 0.87 совершенно не является физическим phiS для набора v=0.2.

И это:

wl
Plot[lam[rr] /. shoot, ...]

тоже неверно: shoot — функция, а не список правил. Именно поэтому Mathematica пишет:

text
{shoot} is neither a list of replacement rules ...

В исправленной версии сначала вычисляется

wl
sol = shoot[phiSsol];

а затем используются непосредственно InterpolatingFunction:

wl
phiIF = sol[[1]];
nuIF  = sol[[3]];
lamIF = sol[[4]];

---

10. Для $e^\nu$ нужна дополнительная нормировка.

Граничное условие горизонта фиксирует

$$e^\nu=a_1(r-r_S)+\ldots,$$

но константа $a_1$ связана с произвольной нормировкой времени. В статье условие (28) задаёт только исчезновение $e^\nu$ и $e^{-\lambda}$ на горизонте. (2207.03944v2.pdf)

Поэтому выбор

wl
nu[r0] == Log[eps]

допустим для интегрирования, но перед сравнением с Fig. 4 надо сделать постоянный сдвиг $\nu$, чтобы метрика была асимптотически Минковской.

Я добавил:

wl
nuShift = -(nuIF[rmax] + lamIF[rmax]);
enu[r_] := Exp[nuIF[r] + nuShift];

---

Сами уравнения (24), (25), (26), (27)

Знаки и коэффициенты в вашей реализации (24) и (25) после решения относительно lambda' и nu' у вас в основном были правильными. Аналогично разложение (26),(27) в rhsVals по $\nu''$ и $\phi''$ математически сделано правильно — там основная ошибка была именно в передаче значений как функций.

И использование (26) здесь не противоречит вашему желанию решать (24),(25),(27): авторы отмечают, что из четырёх уравнений независимы только три из-за тождеств Бианки. (2207.03944v2.pdf) Уравнение (26) можно использовать алгебраически для исключения $\nu''$ из (27).

Исправленный код

Я собрал чистую версию всего основного расчёта: потенциал, (30), (32), правильные (24),(25), исключение вторых производных, shoot, residual, findPhiS и построение трёх графиков. scanResidual туда намеренно не включал.

[Скачать исправленный Wolfram Language код]

Для вашего набора $\mu=0.1,\ell=3,v=0.2$ контрольные значения после исправления должны быть примерно такими:

text
phiS          ~ 6.4*10^-9
phi[50]       ~ 0.2
phi[12]/v     ~ 0.95

Последнее хорошо соответствует фиолетовой кривой Fig. 3: в статье для неё как раз указано $\mu=0.1,\ell=3,v=0.2$. (2207.03944v2.pdf)

Самая важная последовательность проверки сейчас такая: сначала выполните файл целиком, затем отдельно посмотрите phiPrimeS[0.2] — он корректно должен сообщить, что физической вещественной ветви нет; после этого phiSsol = findPhiS[] должен найти малое значение порядка $10^{-9}$, и уже shoot[phiSsol] должен дать InterpolatingFunction, а не $Failed.