# Non-linear numerical solution with variable parameter

I have solved the following system non-linear:

eps = 10^-6; end = 0.2341023826;  g=0.2;
ode1 = a'[r]/r + 4 g (h'[r]^2) (1 - h[r]) == 0;
ode2 = h'[r] (a[r]/r) - (1/4)*(1/g^3) == 0;
bcs = {h[end] == 1, a[eps] == 1};

sol = Map[First[NDSolve[{ode1, ode2, bcs}, {h, a}, r, AccuracyGoal -> 20,
Method -> {"Shooting",
"StartingInitialConditions" -> {h[eps] == # , a[eps] == 1}},
MaxSteps -> 10000]] &, {10^-6}];
Plot[Evaluate[{h[r], h'[r]} /. sol], {r, eps, end}, PlotLegends -> {"h", "h'"}]
Plot[Evaluate[{a[r], a'[r]} /. sol], {r, eps, end}, PlotLegends -> {"a", "a'"}]


This is a physical system and other two condition also must be satisfied; namely, h[0]=0 and a'[end]=0. On the other hand, the value of end should be adjusted according to g. Above it is resolved to g=0.2so that for end=0.2341023826 satisfies the two additional conditions with a good precision, 10^-7 order. However, I could not solve for large values of g. I would like at least until g=2.

Note: In g=1, for example, appear the following error message while trying to adjust the value end:

"NDSolve::ndsz: At r == 1.6338560043295574, step size is effectively zero; singularity or stiff system suspected."

Any help is welcome.

• 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. Aug 16 '18 at 1:05

This problem appears at first glance to have four boundary conditions for a second-order system of ODEs, namely, {h[0] == 0, a[0] == 1, h[end] == 1, a'[end] == 0}. However, ode1 guarantees that a'[end] == 0 when h[end] == 1, reducing the number of independent boundary conditions to three. And, since end can be varied, simply choose it to be the value of r for which h[end] == 1, closing the system. These considerations can be embodied in code as

g = .2;
ode1 = a'[r]/r + 4 g (h'[r]^2) (1 - h[r]) == 0;
ode2 = h'[r] (a[r]/r) - (1/4)*(1/g^3) == 0;
bcs = {h[0] == 0, a[0] == 1};

sol = First@NDSolve[{ode1, ode2, bcs, WhenEvent[h[r] > 1, end = r; "StopIntegration"]},
{h, a}, {r, 0, 5}];
Plot[Evaluate[{h[r], h'[r]} /. sol], {r, 0, end}, PlotLegends -> {"h", "h'"}]
Plot[Evaluate[{a[r], a'[r]} /. sol], {r, 0, end}, PlotLegends -> {"a", "a'"}]


which yields end = 0.234102, and the plots shown in the question. This same method works well for larger values of g, for instance g = 0.4, which yields end = 0.584566 and the plots

Note, however, that a[r] dips toward zero as r approaches end. A small amount of experimentation shows that a[end] == 0 at approximately g = .43135, and both ODEs becomes singular for larger values of g, which accounts for the error message described in the question. Consequently, this nonlinear system has no solutions for larger g`.

• Thanks very much for your answer. Aug 16 '18 at 0:18
• Dear @bbgodfrey, I have a topic already posted for some time link. However, I now need to solve this same problem however, this time, so that a new boundary condition is satisfied. As you have vast experience with solutions of equations differences involving boundary conditions, I would like to know if you could help of any. A colleague from this community has posted an algorithm to help, but has not yet solved the problem. Thanks in advance for your help. Oct 23 '18 at 21:25
• @acsantos I shall see what I can do, but it may not be tonight. Oct 23 '18 at 22:45
• @acsantos The answer by Alex Trounev seems to address well the question that you originally posted, and you accepted it. If you have a new question, please post it separately. Oct 24 '18 at 2:16
• I shall look this evening. Thanks Oct 25 '18 at 12:01