@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] [![enter image description here][1]][1] I hope this helps. [1]: https://i.sstatic.net/1TZr6.jpg