Issue with the NDSolve code

Question

With this procedure, one may determine an eigen-value function $R(a)$ for any given $\Xi$ (say 0, 25, 50, 75, 100)

g[eta1_, eta3_, R_, a_, Xi_] := {f[1], h[1], q[1]} /.     
NDSolve[{D[f[x], {x, 2}] - a^2 f[x] + a (h[x] + R q[x]) == 0,
    D[h[x], { x, 2}] - a^2 h[x] + a Xi x f[x] == 0,
    D[q[x], {x, 2}] - a^2 q[x] + a f[x] == 0,
     f[0] == 0,
f'[0] == 1,
h[0] == eta1,
h'[0] == 0, 
q[0] == 0,
q'[0] == eta3 
}, {f, h, q}, {x, 0, 1}]
s[a_, Xi_] :=FindRoot[g[eta1, eta3, R, a, Xi], {{R, 0}, {eta1, 0}, {eta3, 0}}, 
  MaxIterations -> 100000]

where the condition $f'(0) = 1$ is a normalisation condition. But the problem is that the code doesn't converge always. The output I am looking for is like this

Thanks

Answer 1

Except for simple cases (the only case I can think of is monotonic function at the moment), results of FindRoot are often affected by the starting point for root searching. More starting points, more uncertain, so FindRoot is easier to fail when it's searching for the roots of a set of equations than that of a single equation. The purpose of your code is in fact searching $R$ that makes $f(1)$, $h(1)$, $q(1)$ equal to $0$ for given $a$ and $Ξ$, i.e. $f(x)$, $h(x)$, $q(x)$ will satisfy the boundary condition(B.C.):

$$f(1)=h(1)=q(1)=0$$

Why not use some of these B.C.s instead of the uncertain ones inside NDSolve?

After some trial I found that using $f(1)=q(1)=0$ inside seems to be the best choice, so I modified your code into:

(*
 h[0] == eta1 => f[1] == 0
q'[0] == eta3 => q[1] == 0
*)
g2[R_?NumericQ, a_?NumericQ, Xi_?NumericQ] := 
 h[1] /. NDSolve[{D[f[x], {x, 2}] - a^2 f[x] + a (h[x] + R q[x]) == 0,
                  D[h[x], {x, 2}] - a^2 h[x] + a Xi x f[x] == 0, 
                  D[q[x], {x, 2}] - a^2 q[x] + a f[x] == 0, 
                  f'[0] == 1, h'[0] == 0, q[0] == 0, f[0] == 0, f[1] == 0, q[1] == 0}, 
                 {f, h, q}, {x, 0, 1}]
s2[a_, Xi_] := FindRoot[g2[R, a, Xi], {R, 0}]

The role of those ?NumericQs is to quit the warning NDSolve::ndinnt and ReplaceAll::reps, the output won't be affected even if you take them away.

Let's check the effect of the new definition:

sample = Join[Range[6/10, 2, 1/10], Range[2, 8, 1/4]];
data[Xi_] := {sample, R /. s2[#, Xi] & /@ sample} // Transpose;
rst = data /@ Range[0, 100, 25];
ListLinePlot[rst, PlotRange -> {{0, 8}, {-40, 100}}, Frame -> True, Axes -> False]

Some warnings still generate, but the result matches your figure very well:

And needless to say, much better than your original:

olddata[Xi_] := {sample, R /. s[#, Xi] & /@ sample} // Transpose;
oldrst = olddata /@ Range[0, 100, 25];
ListLinePlot[oldrst, PlotRange -> {{0, 8}, {-40, 100}}, Frame -> True, Axes -> False]