Non-linear numerical solution with variable parameter

Question

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.

Answer 1

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.