# Solving nonlinear 3rd order ODE over range from zero to infinity

Here are the code that I had tried but there is an error

NDSolve[{f'''[y] + f[y] f''[y] + 4 - (f'[y])^2 == 0, f[0] == 0,
f'[0] == 0, f'[y -> ∞] -> 2}, f, {y, 0, 10}]


Here are the question :

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! Nov 28, 2015 at 2:51
• What is the symbol Pr? Note also that NDSolve cannot operate over an infinite domain. So, the upper limit on y must be some finite number. Nov 28, 2015 at 5:00
• I have a feeling this kind of equation has been asked several times in this site. Can you add some background information of the equation so we might be able to solve the problem better? Nov 28, 2015 at 6:39
• @xzczd Certainly. I shall do so tomorrow morning. Just now, I would observe that such problems typically are very difficult to solve numerically, because the solution tracks an unstable attractor of the nonlinear equation. This one is relatively easy, however, because it enters the asymptotic regime so quickly. Nov 28, 2015 at 7:06

f can be computed as follows.

s = NDSolveValue[{f'''[y] + f[y] f''[y] + 4 - (f'[y])^2 == 0, f[0] == 0, f'[0] == 0,
f'[5] == 2}, f, {y, 0, 5}, Method -> "Shooting",
"StartingInitialConditions" -> {f[0] == 0, f'[0] == 0, f''[0] == 3.48}}];
Plot[s[y], {y, 0, 5}, AxesLabel -> {y, f}]
Plot[s'[y], {y, 0, 5}, PlotRange -> {-1, 3}, AxesLabel -> {y, f'}]


Clearly, the solution already is in the asymptotic regime by y = 3.

Computation of θ

With ψ defined as θ', the second ODE can be solved formally as

DSolve[{ψ'[y]/ψ[y] + Pr f[y] == 0, ψ[0] == -1}, ψ[y], y]
{{ψ[y] -> -E^(-Integrate[-(Pr*f[K[1]]), {K[1], 1, 0}] +
Integrate[-(Pr*f[K[1]]), {K[1], 1, y}])}}


To proceed, Pr needs to be specified. For now, set it to unity. Then, the answer just obtained can be written as the following function.

t[y_?NumericQ] := -Exp[- NIntegrate[s[yp], {yp, 0, y},
Method -> {Automatic, "SymbolicProcessing" -> False}]];
Plot[t[y], {y, 0, 5}, AxesLabel -> {y, θ'}]


Finally, θ can be obtained by a second integration, with c5 the constant of integration chosen so that θ[5] == 0, as required by the question.

c5 = NIntegrate[t[yp], {yp, 0, 5}]
(* -1.23953 *)
tt[y_] := NIntegrate[t[yp], {yp, 0, y}] - c5;
Plot[tt[y], {y, 0, 5}, AxesLabel -> {y, θ}, PlotRange -> All]


• How you chose, 'f''[0] == 3.48'?
– zhk
Nov 28, 2015 at 6:38
• @MMM 3.48 is only a guess, so I does not need to be precise. I tried a few guesses until one worked. Nov 28, 2015 at 6:52

@bbgodfrey has explained very well his approach in the answer to this question. But, In my answer, I am trying to solve both the equations simultaneously.

Eqn1 = f'''[x] + f[x] f''[x] + 4 - (f'[x])^2 == 0

Eqn2 = T''[x] + Pr f[x] T'[x] == 0

BC1 = f[0] == 0;

BC2 =  f'[0] == 0;

BC3 = f'[inf1] == 2;

BC4 = T'[0] == -1;

BC5 = T[inf1] == 0;

param1 = {Pr -> 3.97};

inf1 = 5;

Sol1 = NDSolve[{Eqn1, Eqn2, BC1, BC2, BC3, BC4, BC5} /. param1, {f,
T}, {x, 0, inf1}, Method -> {"Shooting",
"StartingInitialConditions" -> {f[0] == 0, f'[0] == 0, f''[0] == 3.48,
T[0] == 0, T'[0] == -10}}]

Plot[{f'[x] /. Sol1, T[x] /. Sol1}, {x, 0, inf1}, PlotRange -> All,
PlotStyle -> {Black, Red}, Frame -> True,
FrameStyle -> Directive[Black, Bold, 12], PlotRange -> All,
Axes -> False]


I hope this helps.