Trouble with the shooting method for boundary value problem of a 4th-order ODE

Question

This is a question about the fluid mechanics equation, which is solved by a similarity solution ($f(t)$, here).

I'm trying to solve the following boundary value problem with shooting method (taken from $(2)(3)(4)$ of this paper):

$f(t)-t f^{\prime}(t)+a\left(f(t)^{3} f^{\prime \prime \prime}(t)\right)^{\prime}=0$

$f(0)=1, f^{\prime}(0)=f^{\prime \prime \prime}(0)=0, f^{\prime \prime}(\infty)=0, f^{\prime}(\infty)=1$

Five boundary conditions are given, in order to determine the unknown parameter $a$.

I choose the ParametricNDSolveValue with the first four boundary conditions, the fifth condition is used conducting the shooting method. Infinity is replaced by t==100000, but there are some errors with the results:

pfun = ParametricNDSolveValue[{f[t] == t f'[t] - a D[f[t]^3 f'''[t], t], 
   f[0] == 1, f'[0] == f'''[0] == 0, f''[100000] == 0}, 
  f'[100000], {t, 0, 100000}, {a}]

FindRoot[pfun[a] - 1, {a, 2}]

Unfortunately, Mathematica gives something like these:

Power::infy: Infinite expression 1/0.^3 encountered.

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

General::stop: Further output of Power::infy will be suppressed during this calculation.

To sum up, my questions are: how can I figure out this ODE to check whether the boundary condition at infinity (in my shooting algorithm I take infinity as t = 100 000) is satisfied? Is my setting wrong? Thanks!

Update:

When I set the xi as x=10, it still doesn't work. There is a singular at t=0. Errors are shown as:

Power::infy: Infinite expression 1/0.^3 encountered.
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.
General::stop: Further output of Power::infy will be suppressed during this calculation.
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.
ParametricNDSolveValue::ndnum: Encountered non-numerical value for a derivative at t$3391 == 0.`.

However, when I change 'a' to '-a', it seems to get a strange answer, which is beyond my expectation. In fact, the value of 'a' should be around 1.22, as stated in an article.

Update2:

The final purpose is to fix this equation:

$f-x f^{\prime}+a \left(f^{R+2}\left|f^{\prime \prime \prime}\right|^{R-1} f^{\prime \prime \prime}\right)^{\prime}=0$ $f(0)=1, f^{\prime}(0)=f^{\prime \prime \prime}(0)=0, f^{\prime \prime}(\infty)=0, f^{\prime}(\infty)=1$

Find 'a' for a specific value of 'R', the prior question is under the condition R=1. I have tried as:

R = 2;
{fsol, asol} = 
  NDSolveValue[{f[t] == 
     t f'[t] - 
      a[t] D[f[t]^(R + 2) (Abs [f'''[t]])^(R - 1)*f'''[t], t], 
    a'[t] == 0, f[0] == 1, f'[0] == f'''[0] == 0, f''[10] == 0, 
    f'[10] == 1}, {f, a}, {t, 0, 10}];
Plot[{fsol[t], asol[t]}, {t, 0, 10}]
y1 = asol[1]

if R=1, y1 = 1.3417, which is corresponding to answer of @xzczd;

When R takes other values, errors appear:

Power::infy: Infinite expression 1/0. encountered.
NDSolveValue::ndnum: Encountered non-numerical value for a derivative at t == 0.`.

So this problem may be difficult to solve，owing to the singular at t==0.

Answer 1

Setting Shooting method outside of *NDSolve* with FindRoot helps, this might be related to the arguable backslide of "Shooting" method.

inf = 10;

{eq, bc} = {f[t] == t f'[t] - a D[f[t]^3 f'''[t], t], {f[0] == 1, f'[0] == 0, 
    f'''[0] == 0, f''[inf] == 0, f'[inf] == 1}};

pfun = ParametricNDSolveValue[{eq, bc[[;; 3]], f''[0] == c}, f, {t, 0, inf}, {a, c}]

parasol = FindRoot[{pfun[a, c]'[inf] == 1, pfun[a, c]''[inf] == 0}, {{a, 2}, {c, 2}}, 
  MaxIterations -> 500]
(* {a -> 1.3417, c -> 0.632144} *)

The result is slightly different from the a = 0.818809^-1 mentioned in the paper, but it's actually a better one, at least when $\infty$ is approximated as 10:

Block[{a = 0.818809^-1}, 
 pfuntst = ParametricNDSolveValue[{eq, bc[[;; 3]], f''[0] == c}, f, {t, 0, inf}, {c}]]

parasoltst = FindRoot[{pfuntst[c]'[inf] == 1}, {{c, 2}}]
(* {c -> 0.661846} *)

{pfuntst[c]'[inf], pfuntst[c]''[inf]} /. parasoltst
(* {1., -0.0041385} *)

{pfun[a, c]'[inf], pfun[a, c]''[inf]} /. parasol
(* {1., 4.93118*10^-15} *)

As you can see, my $f''(\infty)$ is closer to $0$.

The 2 solutions are quite close to each other by the way:

{pfuntst[c] /. parasoltst, pfun[a, c] /. parasol} // ListLinePlot

You may adjust the parameters to e.g. inf = 5 to check the result further.

The code above is tested on v11.3, v12.0.1 and v12.1.1. In v9.0.1 FindRoot spits out a njnum warning and returns unevaluated, which seems to be a bug.