Solution for n == 1
As noted by xzczd, this problem is similar to that in 104170, although more complicated. The answer is easy to obtain for n == 1
, if we happen to know good initial guesses for F'[0]
and G'[0]
. For instance, using {F'[0] == .510, G'[0] == -0.616}
immediately gives

The complete solution is
Plot[Evaluate[{F[η],G[η],H[η]} /. sol], {η, 0, inf}, PlotRange -> All, AxesLabel -> {x, y}]

The question now becomes, how to obtain these good initial guesses. Proceed as follows. First, explicitly set n == 1
in the second and third ODEs, solve for the second derivatives, and Simplify
.
(Flatten@Solve[{F[η]^2 - (G[η] + 1)^2 + (H[η] + ((1 - n)/(1 + n))*η*F[η])*F'[η] ==
(F'[η]^2 + G'[η]^2)^((n - 1)/2)*F''[η] +
F'[η]*((n - 1)*(F'[η]^2 + G'[η]^2)^((n - 3)/2)*(F'[η]*F''[η] + G'[η]*G''[η])),
2*F[η]*(G[η] + 1) + (H[η] + ((1 - n)/(1 + n))*η*F[η])*G'[η] ==
(F'[η]^2 + G'[η]^2)^((n - 1)/2)*G''[η] +
G'[η]*((n - 1)*(F'[η]^2 + G'[η]^2)^((n - 3)/2)*(F'[η]*F''[η] + G'[η]*G''[η]))}]
// Simplify) /. Rule -> Equal
(* {F''[η] == -1 + F[η]^2 - 2 G[η] - G[η]^2 + H[η] F'[η],
G''[η] == 2 F[η] (1 + G[η]) + H[η] G'[η]} *)
Also, the first ODE reduces to
H'[η] == -2*F[η]
Next, derive alternative boundary conditions at large η
. This can be done by linearizing the equations about F[η] == 0
and 1 + G[η] == 0
,
{F''[η] == H[η] F'[η], G''[η] == H[η] G'[η]}
Now, it is evident from the first ODE that H
is approximately constant for large η
, allowing the linearized equations to be integrated once (with constants of integration chosen so that F'[η]
and G'[η]
approach zero at large η
.
{F'[η] == H[η] F[η], G'[η] == H[η] (G[η] + 1)}
The critical advantage of these new boundary conditions at inf
is that they are valid for fairly small values of inf
. Combine these n == 1
ODEs with the new boundary conditions and choose inf == 5
to start.
inf = 5;
sol = NDSolve[{H'[η] == -2*F[η] ,
F''[η] == -1 + F[η]^2 - 2 G[η] - G[η]^2 + H[η] F'[η],
G''[η] == 2 F[η] (1 + G[η]) + H[η] G'[η],
F[0] == 0, G[0] == 0, H[0] == 0,
F'[inf] == H[inf] F[inf], G'[inf] == H[inf] (G[inf] + 1)},
{F, G, H}, {η, 0, inf}, Method -> {"Shooting", "StartingInitialConditions"
-> {F[0] == 0, G[0] == 0, H[0] == 0, F'[0] == .5, G'[0] == -0.7}}];
Flatten@{F'[0], G'[0]} /. sol
(* {0.51028, -0.615914} *)
where the initial guess for {F'[0], G'[0]}
are those in the question. The result, immediately above, is a better set of initial guesses. So, use them with inf == 20
, yielding still better initial guesses (although the ones just above are good enough).
(* {0.510233, -0.615922} *)
which have hardly changed. These are the values on which the initial guesses for the computation at the beginning of the answer are based.
Solutions for various n (corrected)
Solutions for other values of n
can be obtained iteratively by using the initial guesses for n == 1
for a slightly smaller value of n
, and so forth. However, inf
must be somewhat larger to assure that H[η]
is indeed constant near η == inf
Start by defining
inf = 50;
s = ParametricNDSolve[{H'[η] == -2*F[η] - ((1 - n)/(1 + n))*η*F'[η],
F[η]^2 - (G[η] + 1)^2 + (H[η] + ((1 - n)/(1 + n))*η*F[η])*F'[η] ==
(F'[η]^2 + G'[η]^2)^((n - 1)/2)*F''[η] + F'[η]*((n - 1)*(F'[η]^2 + G'[η]^2)^((n - 3)/2)
*(F'[η]*F''[η] + G'[η]*G''[η])),
2*F[η]*(G[η] + 1) + (H[η] + ((1 - n)/(1 + n))*η*F[η])*G'[η] ==
(F'[η]^2 + G'[η]^2)^((n - 1)/2)*G''[η] + G'[η]*((n - 1)*(F'[η]^2 + G'[η]^2)^((n - 3)/2)
*(F'[η]*F''[η] + G'[η]*G''[η])),
F[0] == 0, G[0] == 0, H[0] == 0,
F'[inf] == H[inf] F[inf], G'[inf] == H[inf] (G[inf] + 1)},
{F, G, H, F', G'}, {η, 0, inf}, {n, Fp0, Gp0},
Method -> {"Shooting", "StartingInitialConditions" ->
{F[0] == 0, G[0] == 0, H[0] == 0, F'[0] == Fp0, G'[0] == Gp0}}];
Then, slowly decrease n
, using the computed values of {F'[0], G'[0]}
from a value of n
as initial guesses for the next value of n
. It turns out that step in n
must be quite small as n
approaches 0.6.
guess = {0.510233, -0.615922};
Do[ans = Through[{F, G, H}[n, First@guess, Last@guess]] /. s;
f[n] = First@ans;
guess = D[Through[ans[η]], η][[1 ;; 2]] /. η -> 0;,
{n, 1., .6, -.002}]
Plot[Evaluate@Through[Thread[f[n] /. n -> Range[1., .6, -.1]][η]], {η, 0, 20},
PlotRange -> All, AxesLabel -> {x, y}, AspectRatio -> 5/7]

The total computation requires about 85 sec on my PC.