The problem is to solve the following system ODEs:

I. $\ \ \ \ \ \ \ \dfrac{4}{r}[1+a(r)]\left(\dfrac{dH}{dr}\right)^2+\dfrac{dG}{dr}=0$

II. $\ \ \ \ \dfrac{1}{r}\left(\dfrac{dA}{dr}\right)+f(r)+k^2G(r)=0$

III. $\ \ \ \ \left(\dfrac{dG}{dr}\right)^2+4k^2\left(G(r)+\dfrac{f(r)}{k^2}\right)^2\left(\dfrac{dH}{dr}\right)^2-(1.6)H(r)\left(\dfrac{dH}{dr}\right)^2=0$

IV. $\ \ \ \ \dfrac{d^2f}{dr^2}+\dfrac{1}{r}\dfrac{df}{dr}+\dfrac{1}{r}\dfrac{dA}{dr}-4k^2f(r)\left(\dfrac{dH}{dr}\right)^2=0$

where $k>0$. The boundary conditions are

$H(0)=1, \ \ A(0)=0, \ \ G(L)=0, \ \ f'(0)=0, \ \ f(L)=0 \ \ [\text{and} \ \ f'(L)=0]$

where $L$ is the boundary, which must be adjusted in order to satisfy the boundary conditions above for each value of parameter $k$ chosen. Note that we have 6 (six) boundary conditions which must be satisfied, but we need only 5 (five). This sixth condition needs to be guaranteed in some way and this is my difficulty. The condition $f'(L)=0$ is necessary due to be a physical problem. Further, this physical problem requires also:

$H(L)=H'(L)=0, \ \ A'(L)=0$.

Thus, assuming these conditions mentioned, a brief analysis of the ODEs allows conclude that:

The ODE I ensures $H'(0)=0$ and $G'(L)=0$ wiht 1 (one) free parameter $a(L)=const.$

The ODE II ensures $A'(0)=0$ and 2 (two) free parameters $G(0)=const.$ and $f(0)=const.$

The ODE III ensures $G'(0)=0$.

It was possible to find a numerical solution using the shooting method as follow, but already advance that in this case the boundary condition A'[L]=0 is not satisfied as I want. It follow that,

ode1 = 4 (h'[r] ^2) (a[r] + 1)/r == -G'[r];
ode2 = a'[r]/r == -(k^2) (G[r]) - f[r];
ode3 = (G'[r])^2 + 4 (k^2) (h'[r]^2) (G[r] + f[r]/k^2)^2 == (8/5) h[r] (h'[r]^2);
ode4 = f''[r] + f'[r]/r + (a'[r]/r) == 4 (k^2) (f[r]) ((h'[r])^2);

Simplifying ode1 and ode3 can avoid problems with complex numbers. This is done using Solve[{ode1, ode3}, {h'[r], G'[r]}] // FullSimplify. The new set of equations is reduced to form

ode13 = {Derivative[1][h][r] == -((r Sqrt[2k^2 h[r] - 5 (f[r] + k^2 G[r])^2])/(2 Sqrt[5]k(1+a[r]))),Derivative[1][G][r] == (r(-2 k^2 h[r] + 5 ( f[r] + k^2 W[r])^2))/(5k^2 (1 + a[r]))};
ode2 = a'[r]/r == -(k^2) (G[r]) - f[r];
ode4 = f''[r] + f'[r]/r + (a'[r]/r) == 4 (k^2) (f[r]) ((h'[r])^2);

Applying the algorithm for k=1, as example, we have:

eps = 10^-10; L = 2.667; k = 1;
ode13 = {Derivative[1][h][r] == -((r Sqrt[2k^2 h[r] - 5 (f[r] + k^2 G[r])^2])/(2 Sqrt[5]k(1+a[r]))),
Derivative[1][G][r] == (r(-2 k^2 h[r] + 5 ( f[r] + k^2 G[r])^2))/(5k^2 (1 + a[r]))};
ode2 = a'[r]/r == -(k^2) (G[r]) - f[r];
ode4 = f''[r] + f'[r]/r + (a'[r]/r) == 4 (k^2) (f[r]) ((h'[r])^2);
bcs = {h[eps] == 1, a[eps] == 0, G[L] == 0, f'[eps] == 0, f[L] == 0};
sol = NDSolve[{ode13, ode2, ode4, bcs}, {a, h, f, G}, {r, eps, L},
Method -> {"Shooting","StartingInitialConditions" -> {f[0] == -0.340181, G[0] == 0.736905}}];
Plot[Evaluate[Re[{h[r], h'[r]}] /. sol], {r, eps, L}, PlotLegends -> {"h", "h'"}, PlotRange -> All]
Plot[Evaluate[Re[{a[r], a'[r]}] /. sol], {r, eps, L}, PlotLegends -> {"a", "a'"}, PlotRange -> All]
Plot[Evaluate[Re[{f[r], f'[r]}] /. sol], {r, eps, L}, PlotLegends -> {"f", "f'"}, PlotRange -> All]
Plot[Evaluate[Re[{G[r], G'[r]}] /. sol], {r, eps, L}, PlotLegends -> {"G", "G'"}, PlotRange -> All]

enter image description here

In this case, the values of f[0] and G[0] were adjusted together with the value of L in order to satisfy the boundary conditions. With the exception of f'[L]=0, all the other boundary conditions has been satisfied.

Then, main question is:

Is it possible to obtain a numerical solution by assuring the boundary condition f'[L]=0 in some way together with the others already reached in the solution shown above?

Remark: Note that in the above solution, it was only the condition f[L]=0 was used and not f'[L]=0. I tried to substitute one for the other but not works, that is, it does not satisfy the problem.

  • 2
    $\begingroup$ It is necessary to give a link to the original topic and all solutions mathematica.stackexchange.com/questions/180269/… $\endgroup$ Oct 24, 2018 at 6:45
  • 1
    $\begingroup$ For each additional boundary condition, a parameter is required. For example, you can pick k for which f'[L]=0 $\endgroup$ Oct 24, 2018 at 6:57
  • 1
    $\begingroup$ How do you obtain h[L] == 0 mathematically from the equations? $\endgroup$
    – bbgodfrey
    Oct 27, 2018 at 1:19
  • 1
    $\begingroup$ @acsantos With Normal@Series[eq, {r, r0, 0}] /. {f[r0] -> 0, G[r0] -> 0} with r0 an arbitrary value of r > 0 yields a'[r0] -> 0 but not h[r0] -> 0, unless r0 is equal to the particular value of L given in your question. In other words, you appear to be choosing L to cause h[L] to vanish. However, there appears to be a branch point in the solution at that value of L, and integrating past is does not seem possible. In any case, I have spent some time attempting to find a way to obtain f'r0] -> 0` but without success. Probably, it cannot be done. $\endgroup$
    – bbgodfrey
    Oct 28, 2018 at 2:36
  • 1
    $\begingroup$ @acsantos By the way, Solve[{ode1, ode3}, {h'[r], G'[r]}] has three solutions. Are you sure that you have chosen the correct one? $\endgroup$
    – bbgodfrey
    Oct 28, 2018 at 2:38

1 Answer 1


In principle, it should be possible to use bisection search to find a suitable parameter L. Define

ClearAll[F, L];
F[L_?NumericQ] := 
 f'[L] /. NDSolve[{ode13, ode2, ode4, bcs}, {a, h, f, G}, {r, eps, L},
     Method -> {"Shooting", 
      "StartingInitialConditions" -> {f[0] == -0.340181, 
        G[0] == 0.736905}}][[1]]

and run

FindRoot[F[L], {L, L0, L1}, Method -> "Brent"]

where L0 and L1 have to bracket a solution, i.e., F[L0] and F[L1] need to have opposite signs. But in your particular case, this seems to be impossible: NDSolve produces bogus for values greater than L = 2.667 and for values below 2.667, we have

LList = Subdivide[0.0000001, 2.667, 100];
data = F /@ LList;
ListLinePlot[Transpose[{LList, data}], AxesLabel -> {"L", "f'[L]"}, 
 PlotRange -> {0, Max[data] 1.1}]

enter image description here

This sort of implies that the only solution to your problem is L=0.

  • $\begingroup$ Dear @HenrikShumacher, Could I take the derivative of any of the equations, for example Eq. II, in order to introduce the boundary condition f'[L]=0 and, thus, try to solve the problem? Is it acceptable to do this? $\endgroup$
    – Andrey
    Oct 25, 2018 at 1:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.