How to solve this nonlinear differential equation?

Question

The given nonlinear differential equation is

y'''[t]+(y[t]*y''[t])+y[t]'^2-1=0

with boundary conditions {y[0]=0,y'[0]=0 and y'[t]->1 as t->Infinity.

Answer 1

Answer substantially revised.

Numerical Solution

This question looks similar to 100659, so one might expect to solve it in the same way. For instance, write

tmax = 1000; 
s = ParametricNDSolveValue[{y'''[t] + y[t]*y''[t] + y'[t]^2 - 1 == 0, 
    y[0] == 0, y'[0] == 0, y''[0] == c}, y, {t, 0, tmax}, {c}];

and then apply FindRoot to s[c][tmax] == 1 to find the value of c for which y'[tmax] == 1. NDSolve with the option Method -> {Shooting, ...} would be expected to have the same result. However, for this question a range of c satisfy the outer boundary condition. This can be seen by plotting solutions for several values of c.

Plot[{s[-1.08637576][t], s[-1][t], s[0][t], s[1][t], s[3][t], s[10][t]}, 
    {t, 0, tmax/50}, AxesLabel -> {t, y}, LabelStyle -> Directive[Black, Bold, Medium]] 
Plot[{s[-1.08637576]'[t], s[-1]'[t], s[0]'[t], s[1]'[t], s[3]'[t], s[10]'[t]}, 
    {t, 0, tmax/50}, AxesLabel -> {t, y'}, LabelStyle -> Directive[Black, Bold, Medium], 
    PlotRange -> {-1, 5}]
{s[-1.08637576]'[tmax], s[-1]'[tmax], s[0]'[tmax], s[1]'[tmax], s[3]'[tmax], s[10]'[tmax]}

(* {1., 1., 1., 1., 1.00001, 1.00005} *)

For values of c smaller than about -1.08637576, the solution is singular at finite t.

Symbolic Solution

Interestingly, this question also can be solved symbolically. First, observe that the ODE can be written as

D[y'[t] + y[t]^2/2, {t, 2}] -1 == 0

Obviously, this can be twice integrated. Set z[t] == y'[t] + y[t]^2/2, and

(Flatten@DSolve[{z''[t] == 1, z[0] == 0}, z[t], t] // Apart) /. C[2] -> c
(* {z[t] -> c t + t^2/2} *)

So, the ODE can be rewritten as a first-order differential equation,

y'[t] + y[t]^2/2 == t^2/2 + c t

with the remaining boundary condition, y[0] == 0. Then, DSolve yields

Simplify@DSolveValue[{y'[t] + y[t]^2/2 == t^2/2 + c t, y[0] == 0}, y[t], t]

(-(c + t) ParabolicCylinderD[1/4 (-2 - c^2), I (c + t)] (c ParabolicCylinderD[1/4 (-2 + c^2), c] - 2 ParabolicCylinderD[1/4 (2 + c^2), c]) - 2 I ParabolicCylinderD[1/4 (2 - c^2), I (c + t)] (c ParabolicCylinderD[1/4 (-2 + c^2), c] - 2 ParabolicCylinderD[1/4 (2 + c^2), c]) + (c ParabolicCylinderD[ 1/4 (-2 - c^2), I c] + 2 I ParabolicCylinderD[1/2 - c^2/4, I c]) ((c + t) ParabolicCylinderD[1/4 (-2 + c^2), c + t] - 2 ParabolicCylinderD[1/4 (2 + c^2), c + t]))/((c ParabolicCylinderD[1/4 (-2 - c^2), I c] + 2 I ParabolicCylinderD[1/2 - c^2/4, I c]) ParabolicCylinderD[ 1/4 (-2 + c^2), c + t] + ParabolicCylinderD[1/4 (-2 - c^2), I (c + t)] (c ParabolicCylinderD[1/4 (-2 + c^2), c] - 2 ParabolicCylinderD[1/4 (2 + c^2), c]))

Plotting this expression yields the same curves as above with the denominator becoming singular for c about equal to -1.08637574066271760496838769091.