# Issue with the NDSolve code

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

• Please add some code to demonstrate where your implementation allegedly fails. Oct 28, 2013 at 13:57
• Well, for s[7.8,100], R=28 and for s[8,100], R=22, which is not what we can see from the above fig.
– zhk
Oct 28, 2013 at 14:26
• Please update your question accordingly and give us a simple template to check against the figure. Oct 28, 2013 at 14:46
• I noticed my mistake and rewrited my answer, I think now it's just what you want. Have a look! Nov 4, 2013 at 13:11

## 1 Answer

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]


• lovely! thx mate
– zhk
Nov 4, 2013 at 17:05