# Solving nonlinear ODE with boundary conditions

I am trying to solve the below fucntions using a NDsolve, but I can't get the the converges of them at the certain initial variables.I'm using Mathematica 11.3. I appreciate any help.

Z=800;
g= 0.023800000000000000000;
k2= 0.000194519;
R= 1.5472;
ytest0= -13.911917694213733;
ϵ = \$MachineEpsilon ;
y1[ytest_?NumericQ] :=
NDSolve[
{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]] , y[1] == ytest,
y'[ϵ] == 0}, y, {r, ϵ, 1},
Method -> {"StiffnessSwitching", "NonstiffTest" -> False}];
y2[ytest_?NumericQ] :=
NDSolve[
{y''[r] + 2 y'[r]/r == κ2 Sinh[y[r]],
y[1] == ytest, y'[R] == 0}, y, {r, 1, R},
Method -> {"StiffnessSwitching", "NonstiffTest" -> False}];
y1Try[ytest_?NumericQ] := y'[1] /. y1[ytest];
y2Try[ytest_?NumericQ] := y'[1] /. y2[ytest];
f = ytest /. FindRoot[y1Try[ytest] - y2Try[ytest]==-Zg, {ytest, ytest0}]

• @zhk In the way that every ODE is also a PDE... ;) – Henrik Schumacher Sep 16 '18 at 12:35
• Cross Post – Mariusz Iwaniuk Sep 16 '18 at 12:57
• The discontinuity Z/g in y'[r] at r == 1 is enormous. Are you sure that this value is correct? – bbgodfrey Sep 16 '18 at 18:19
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Sep 17 '18 at 1:17

## 1 Answer

The quantity κ2 in the question is undefined. However, the cross-posting identified by Mariusz Iwaniuk in a comment above suggests that k2 actually is meant. With this change, the computation is equivalent to solving the ODE, {y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], y'[R] == 0, y'[ϵ] == 0} with a discontinuity in y'[r] equal to Z g at r == 1. With this modification to the ODE, it can be recast as

Z = 800; g = Rationalize[0.0238, 0];
k2 = Rationalize[ 0.000194519, 0];
ϵ = 10^-4; R = Rationalize[1.5472, 0];
ps = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]],
y[ϵ] == y0, y'[ϵ] == 0, WhenEvent[r == 1, y'[r] -> y'[r] + Z g]}, {y, y'},
{r, ϵ, R}, {y0}, Method -> "StiffnessSwitching", WorkingPrecision -> 30];
sol = FindRoot[Last[ps[y0]][R], {y0, -11.25, -11.}, Evaluated -> False][[1, 2]]
(* -11.1979 *)

Plot[First[ps[sol]][r], {r, ϵ, R}, AxesLabel -> {r, y}, ImageSize -> Large,
LabelStyle -> {Black, Bold, Medium}]


Plot[Last[ps[sol]][r], {r, ϵ, R}, PlotRange -> All, AxesLabel -> {r, y'},
ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]


The computation takes only seconds.

Addendum

In response to a comment below, the solution for k2 = Rationalize[0.05716690559951900467, 0]; is

sol = Quiet@FindRoot[Last[ps[y0]][R], {y0, -12, 0}, Evaluated -> False][[1, 2]]
(* -5.51467 *)


The corresponding plot of y looks very much like that above but shifted approximately by the shift in y0.

The corresponding plot of y' also looks very much like that above. Note that the FindRoot search range has been expanded in this computation, and Quiet has been prepended to avoid printing the numerous errors that FindRoot encounters before finding the correct answer.

Second Addendum

Varying k2, as requested by the OP in a comment below, is accomplished as follows.

psk = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == k Sinh[y[r]],
y[ϵ] == y0, y'[ϵ] == 0, WhenEvent[r == 1, y'[r] -> y'[r] + Z g]}, {y, y'},
{r, ϵ, R}, {y0, k}, Method -> "StiffnessSwitching", WorkingPrecision -> 30];
ytablog = Table[{Exp[k], Quiet@FindRoot[Last[psk[y0, Exp[k]]][R], {y0, -12, 0},
Evaluated -> False][[1, 2]]}, {k, -9, 0, 1/10}];
ListLogLinearPlot[ytablog, PlotRange -> All, Joined -> True, AxesLabel -> {k2, y[0]},
ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]


A rough fit to this curve is Log[k2] - 2.655. Varying R can be accomplished in a similar manner.

Interestingly, varying R between 5/4 and 2 has minimal effect on the plot just above. Below are plots of y[r] for R = Range[5/4,2,1/4] - {"Blue", "Orange", "Green", "Red"}.

• Thank you bbgodfrey, There is a mistake in my code the discontinuity (Zg) in y'[r] at r == 1 not (Z/g) sorry for that. From that, I am looking for discontinuity in region [20 to 30], when using your way for 23 given me a wrong result even though any values from 5 to 25 gave a wrong result.So, I am not sure that (28) number comes from where? it should be based on the initial values. Does it come from multiply (Z by g)? – aluuzz Sep 17 '18 at 4:51
• @user60211 I used 28, because it yielded a value of y[1] approximately equal to your ytest0. Rerunning the code with 19.04 is easy, but it is very late here. I shall do so tomorrow morning. By the way, what do you mean by "wrong result" in your comment? – bbgodfrey Sep 17 '18 at 5:00
• I meant the graph is not reasonable. I am looking for something similar what you got for 28. Ok I am waiting you and thanks for your help. – aluuzz Sep 17 '18 at 5:08
• @user60211 Redoing the computation and modifying my answer for Z g instad of 28` did not take long. I hope it meets your need. – bbgodfrey Sep 17 '18 at 5:36
• very nice! I got what I want. Thank you so much. Last question if I increase the k2 value let say k2 = 0.05716690559951900467 my code is not running. Can you help me in that. By the way, Can I plot the data from 0 to 1? – aluuzz Sep 17 '18 at 6:43