Kernel crashes after long iterations

Question

I am trying to get a numerical root of a function defined as below. But after some large recursion the kernel get crashed. But for relatively low recursion(10000) that does not make any problem. But I need to do the recursion for about 200000 time. Here I am attaching the code. Anyone please help me to overcome this problem.

RootFinderDo[(\[Omega]_)?NumberQ, L_, MAX_, r_, M_, Q_, dQdr_] := (ClearAll[A, B, G, \[Alpha], \[Beta], \[Gamma], \[Delta], S]; \[Alpha][n_] := ((1 - 2*M)/r)*(n*(n + 1)); \[Beta][n_] := -2*n*(I*(\[Omega]/M)*r + (1 - (3*M)/r)*n); \[Gamma][n_] := (1 - (6*M)/r)*n*(n - 1) + (6*M)/r - L*(L + 1); \[Delta][n_] := ((2*M)/r)*(n - 3)*(n + 1);A[0] = -1; B[0] = r*(dQdr/Q + (I*(\[Omega]/M)*r)/(r - 2*M)); A[1] = \[Alpha][1]; B[1] = \[Beta][1]; G[1] = \[Gamma][1]; A[n_] := A[n] = \[Alpha][n]; G[n_] := G[n] = \[Gamma][n] - (B[n - 1]*\[Delta][n])/G[n - 1]; B[n_] := B[n] = \[Beta][n] - (A[n - 1]*\[Delta][n])/G[n - 1]; S[MAX] = G[MAX]/B[MAX];Do[S[n] = G[n + 1]/(B[n + 1] - A[n + 1]*S[n + 1]), {n, MAX - 1, 0, -1}]; RR = B[0] - A[0]*S[0])

\[Omega]guess = 0.0065

Block[{$RecursionLimit = 20000}, FindRoot[RootFinderDo[\[Omega], 2, 15000, 1.2579481679118418*^6, 196567.4968628099, -761321.9343226571 - 24959.124711629214*I, -1.9220363777903136 - 0.06301190533310615*I] == 0, {\[Omega], \[Omega]guess}]]

Answer 1

There's no need for recursive definitions, you can do everything iteratively:

RootFinderDo[ω_?NumberQ, L_, MAX_, r_, M_, Q_, dQdr_] :=
  Module[{A, B, G, α, β, γ, δ, S},
    α[n_] = ((1 - 2*M)/r)*(n*(n + 1));
    β[n_] = -2*n*(I*(ω/M)*r + (1 - (3*M)/r)*n);
    γ[n_] = (1 - (6*M)/r)*n*(n - 1) + (6*M)/r - L*(L + 1);
    δ[n_] = ((2*M)/r)*(n - 3)*(n + 1);
    A[0] = -1;
    A[n_] = α[n]; 
    B[0] = r*(dQdr/Q + (I*(ω/M)*r)/(r - 2*M));
    B[1] = β[1];
    G[1] = γ[1];
    Do[B[n] = β[n] - (A[n - 1]*δ[n])/G[n - 1];
       G[n] = γ[n] - (B[n - 1]*δ[n])/G[n - 1], {n, 2, MAX}];
    S[MAX] = G[MAX]/B[MAX]; 
    Do[S[n] = G[n + 1]/(B[n + 1] - A[n + 1]*S[n + 1]), {n, MAX - 1, 0, -1}];
    B[0] - A[0]*S[0]]

Notice that I used = everywhere, not :=, to make the definitions stick immediately.

Also notice the use of Module, which localizes the variables and makes it unnecessary to clean them manually or to take care of the $HistoryLength manually.

On the "return trip" there is actually no need to store the values of S[n]; we can just keep the "latest" value of S in memory:

RootFinderDo[ω_?NumberQ, L_, MAX_, r_, M_, Q_, dQdr_] :=
  Module[{A, B, G, α, β, γ, δ, S},
    α[n_] = ((1 - 2*M)/r)*(n*(n + 1));
    β[n_] = -2*n*(I*(ω/M)*r + (1 - (3*M)/r)*n);
    γ[n_] = (1 - (6*M)/r)*n*(n - 1) + (6*M)/r - L*(L + 1);
    δ[n_] = ((2*M)/r)*(n - 3)*(n + 1);
    A[0] = -1;
    A[n_] = α[n]; 
    B[0] = r*(dQdr/Q + (I*(ω/M)*r)/(r - 2*M));
    B[1] = β[1];
    G[1] = γ[1];
    Do[{B[n], G[n]} = {β[n], γ[n]} - {A[n - 1], B[n - 1]}*δ[n]/G[n - 1],
       {n, 2, MAX}];
    S = G[MAX]/B[MAX]; 
    Do[S = G[n + 1]/(B[n + 1] - A[n + 1]*S), {n, MAX - 1, 0, -1}];
    B[0] - A[0]*S]

This saves about 10% of time.

I don't think there's a root where you expect it:

ReImPlot[RootFinderDo[ω, 2, 15000, 1.2579481679118418*^6, 196567.4968628099,
                      -761321.9343226571 - 24959.124711629214*I,
                      -1.9220363777903136 - 0.06301190533310615*I],
         {ω, 0, 0.01}]