Need help in getting a numerical solution to a nonlinear ODE

I am trying to solve a nonlinear ODE BY applying a StiffnessSwitching method, but when I try to find the root of my equation it gives me a error message. It seems to me the problem is that Mathematica 11 passes a variable's quantity instead of passing a number in the FindRoot.

Does anyone know how to fix this problem?

Subscript[\[Lambda], B] = 0.71410720;
a = 25;
Z = 1000;
Subscript[c, 0] = .2;
\[Eta] = 0.3;
\[Gamma] = Subscript[\[Lambda], B]/a;
R = \[Eta]^(-1/3);
Subscript[n, m] = 3 \[Eta]/(4 \[Pi]);
Subscript[n, 0] = Subscript[c, 0] Avogadro a^3;
\[Kappa]2 = 8 \[Pi] \[Gamma] Subscript[n, 0] ;
\[Kappa] = Sqrt[\[Kappa]2] ;


ψtest is 5

 in[ψtest_?NumericQ] :=
NDSolve[
{ψ''[r] + 2 ψ'[r]/r == κ2 Sinh[ψ[r]] + 3 Z γ, ψ[1] == ψtest,
ψ'[ϵ] == 0}, ψ, {r, ϵ, 1},
Method -> {"StiffnessSwitching", "NonstiffTest" -> False}];
ψout[ψtest_?NumericQ] :=
NDSolve[
{ψ''[r] + 2 ψ'[r]/r == κ2 Sinh[ψ[r]],
ψ[1] == ψtest, ψ'[R] == 0}, ψ, {r, 1, R},
Method -> {"StiffnessSwitching", "NonstiffTest" -> False}];
ψinTry[ψtest_?NumericQ] := ψ'[1] /. ψin[ψtest];
ψoutTry[ψtest_?NumericQ] := ψ'[1] /. ψout[ψtest];
a = ψtest /. FindRoot[ψinTry[ψtest] == ψoutTry[ψtest], {ψtest, ψtest0}]

• I don't see any values given for Kappa2, Z or Gamma. – Bill Watts Dec 6 '17 at 0:26
• These are just constants i did not include them in my post – user53936 Dec 6 '17 at 9:22
• We can't help debug without knowing those values. – Bill Watts Dec 6 '17 at 18:42
• I updated my code – user53936 Dec 6 '17 at 21:43
• Still missing values for \[Psi]test0, [Epsilon], maybe others, and I assume you mean \[Psi]in instead of in. – Bill Watts Dec 6 '17 at 23:27

I think with the date you are using, your FindRoot will always fail. I see you are trying to align the slopes of wave functions that are refusing to align. I broke your code into pieces to get a better look:

\[Epsilon] = 2/10;
\[Psi]test = -2;

NDSolve[{\[Psi]''[r] + 2 \[Psi]'[r]/r == \[Kappa]2 Sinh[\[Psi][r]] + 3 Z \[Gamma],
\[Psi][1] == \[Psi]test, \[Psi]'[\[Epsilon]] == 0},
\[Psi], {r, \[Epsilon], 1},
Method -> {"StiffnessSwitching", "NonstiffTest" -> False},
WorkingPrecision -> 30];

\[Psi]1[r_] = \[Psi][r] /. %[[1]];

NDSolve[{\[Psi]''[r] + 2 \[Psi]'[r]/r == \[Kappa]2 Sinh[\[Psi][r]]
\[Psi][1] == \[Psi]test,
\[Psi]'[R] == 0}, \[Psi], {r, 1, R},
Method -> {"StiffnessSwitching", "NonstiffTest" -> False},
WorkingPrecision -> 30];

\[Psi]2[r_] = \[Psi][r] /. %[[1]];
p1 = Plot[\[Psi]1[r], {r, \[Epsilon], 1}];
p2 = Plot[\[Psi]2[r], {r, 1, R}];

Show[p1, p2, PlotRange -> All]


You can play around with the numbers, but with the given data, I don't think there is any \[Psi]test value that will equate the 2 slopes at psi = 1, and as epsilon gets closer to zero, the first wavefunction NDSolve gets harder to find any solution.

• Ok , but my code is working very well with Mathematica version 9. what did you plot is not the write solution of my problem.I am wondering is there different between the Mathematica 11 or 9 in NDSolve ? – user53936 Dec 7 '17 at 17:25
• OK if M9 works, that's great. Don't have 9 but M10 does the same thing. – Bill Watts Dec 7 '17 at 20:42
• Yeah, But what is the problem why it did not work in M11? I have to run it in M11 – user53936 Dec 7 '17 at 21:10