Using NDSolve to solve three coupled ODEs

This seventh-order system of ODEs is highly nonlinear but nonetheless can be solved iteratively as follows. For simplicity, name the set of three equations eq, the set of four η = 0 boundary conditions bc0, and the set of three η = ub boundary conditions bc1. (ub is the upper bound of the computation, given by 50 in the question.) Then,

soln = NDSolve[{eq, bc0, bc1}, {f, θ, ϕ}, {η, 0, ub}] 

which fails with the error,

NDSolve::ndsz: At η == 1.51819759348527`, step size is effectively zero; singularity or stiff system suspected.

Using Method -> "StiffnessSwitching" does not help.

Note that the boundary condition f'[ub] == 1 suggests that f[ub] may be quite large compared to boundary values at η = 0. The substitution

{eqg, bc0g, bc1g} = Simplify[{eq, bc0, bc1} /. f -> Function[{η}, g[η] + η]] 

eliminates this potential issue but, unfortunately, does not avoid the error described above. I next expanded eqg for large η, which indicates that eqg[[1]] is approximately independent of {θ,ϕ} there. This suggests formally decoupling the first equation, which then can be solved to obtain,

gg[1] = NDSolveValue[{eqg[[1]] /. θ -> Function[η, 0], bc0g[[1 ;; 2]], bc1g[[1]]}, g[η], {η], 0, ub}] // Flatten; Plot[gg[1], {η, 0, ub}, PlotRange -> All, AxesLabel -> {η, g}, LabelStyle -> {10, Bold, Black}] 

which indicates that ub = 20 is plenty large, at least for this level of approximation. Corresponding values of {θ,ϕ} are obtained from

θϕ[1] = NDSolveValue[{eqg[[2 ;; 3]] /. g -> gg[1][[0]], bc0g[[3 ;; 4]], bc1g[[2 ;; 3]]}, {θ[η], ϕ[η]}, {η, 0, ub}] // Flatten; 

These values then can be used to obtain a better approximation for g[2], and so forth. After three iterations,

Plot[{gg[1], gg[2], gg[3]}, {η, 0, ub}, PlotRange -> All, AxesLabel -> {η, g}, LabelStyle -> {10, Bold, Black}] 

gg[2] and gg[3] are indistinguishable to the eye in this plot and differ numerically by about one part in 1000 for large η. The complete solution after the third iteration is

Plot[{gg[3], θϕ[3]}, {η, 0, ub}, PlotRange -> All, AxesLabel -> {η, "g,θ,ϕ"}, LabelStyle -> {10, Bold, Black}] 

Not an answer, but too long for a comment.

The approximation at 50 is giving you the stiffness/singularity problem.

Trying a different upper bound:

ub = 1.5; soln = NDSolve[{Exp[-\[Alpha] \[Theta][\[Eta]]]*(1 + We^2 (f''[\[Eta]])^2)^((n - 3)/ 2)*(-\[Alpha] \[Theta]'[\[Eta]] f''[\[Eta]] (1 + We^2 (f''[\[Eta]])^2) + (1 + n We^2 (f''[\[Eta]])^2) f'''[\[Eta]]) + (b/2 \[Eta] + f[\[Eta]]) f''[\[Eta]] + (b - M - f'[\[Eta]]) f'[\[Eta]] + (M - b + 1) == 0, (1/Pr) (\[CurlyEpsilon] \[Theta]'[\[Eta]]^2 + (1 + \ \[CurlyEpsilon] \[Theta][\[Eta]]) \[Theta]''[\[Eta]]) + (f[\[Eta]\ ] \[Theta]'[\[Eta]] - f'[\[Eta]] \[Theta][\[Eta]]) - Ec (f''[\[Eta]])^2 Exp[-\[Alpha] \[Theta][\[Eta]]]*(1 + We^2 (f''[\[Eta]])^2)^((n - 1)/2) + (b/ 2) (\[Theta][\[Eta]] - \[Eta] \[Theta]'[\[Eta]]) - \ \[CapitalLambda]1 \[Theta]'[\[Eta]]^2 - \[CapitalLambda]2 \[Theta]'[\ \[Eta]] \[Phi]'[\[Eta]] == 0, \[Phi]''[\[Eta]] + Pr Le (f[\[Eta]] \[Phi]'[\[Eta]] - f'[\[Eta]] \[Phi][\[Eta]] - b/2 (\[Phi][\[Eta]] - \[Eta] \[Phi]'[\[Eta]])) + \ (\[CapitalLambda]1/\[CapitalLambda]2) \[Theta]''[\[Eta]] == 0, f[0] == -f0, f'[0] == 0, \[Theta][0] == 1, \[Phi][0] == 1, f'[ub] == 1, \[Theta][ub] == 0, \[Phi][ub] == 0}, {f, \[Theta], \[Phi]}, {\[Eta], 0, ub}]; plot1 = Plot[Evaluate[{f[\[Eta]]} /. soln], {\[Eta], 0, ub}, Frame -> True, PlotLegends -> {"\[Phi](\[Eta])"}, PlotLabel -> "Concentration Profile \[Phi](\[Eta])"] 

We can use NestList to find a solution.

We start iteration with \[Phi]=Function[\[Eta], Exp[-\[Eta] ]], \[Theta]=Function[\[Eta], Exp[-\[Eta] ]](ansatz fullfills boundary conditions) and solve the first ode for f using a selfmade shootingmethod which enforces f'[\[Eta]max]==1.

Knowing iteration for f we can solve second&third ode.

\[Eta]max = 10; (* sufficient large *) solUN = NestList[ (solf = ParametricNDSolveValue[{(*First Equation:Momentum*) Exp[-\[Alpha] #[[2]][\[Eta]]]*(1 + We^2 (f''[\[Eta]])^2)^((n - 3)/ 2)*(-\[Alpha] #[[2]]'[\[Eta]] f''[\[Eta]] (1 + We^2 (f''[\[Eta]])^2) + (1 + n We^2 (f''[\[Eta]])^2) f'''[\[Eta]]) + (b/ 2 \[Eta] + f[\[Eta]]) f''[\[Eta]] + (b - M - f'[\[Eta]]) f'[\[Eta]] + (M - b + 1) == 0 , f[0] == -f0, f'[0] == 0, f''[0] == fss0 (*Boundary Conditions as \[Eta]\[RightArrow]\[Infinity], approximated at \[Eta]=50*)(*f'[\[Eta]max]\[Equal]1*) }, f, {\[Eta], 0, \[Eta]max}, {fss0}, Method -> "StiffnessSwitching"] ; valfss0 = fss0 /. NMinimize[(solf[fss0]'[\[Eta]max] - 1)^2, fss0][[2]] // Quiet ; ff = solf[valfss0]; (* iterierte Lösung f*) sol\[Theta]\[Phi] = NDSolveValue[{ (*Second Equation: Energy*)(1/ Pr) (\[CurlyEpsilon] \[Theta]'[\[Eta]]^2 + (1 + \ \[CurlyEpsilon] \[Theta][\[Eta]]) \[Theta]''[\[Eta]]) + (ff[\[Eta]] \ \[Theta]'[\[Eta]] - ff'[\[Eta]] \[Theta][\[Eta]]) - Ec (ff''[\[Eta]])^2 Exp[-\[Alpha] \[Theta][\[Eta]]]*(1 + We^2 (ff''[\[Eta]])^2)^((n - 1)/2) + (b/ 2) (\[Theta][\[Eta]] - \[Eta] \[Theta]'[\[Eta]]) - \ \[CapitalLambda]1 \[Theta]'[\[Eta]]^2 - \[CapitalLambda]2 \[Theta]'[\ \[Eta]] \[Phi]'[\[Eta]] == 0,(*Third Equation: Concentration*)\[Phi]''[\[Eta]] + Pr Le (ff[\[Eta]] \[Phi]'[\[Eta]] - ff'[\[Eta]] \[Phi][\[Eta]] - b/2 (\[Phi][\[Eta]] - \[Eta] \[Phi]'[\[Eta]])) + (\ \[CapitalLambda]1/\[CapitalLambda]2) \[Theta]''[\[Eta]] == 0,(*Boundary Conditions at \[Eta]= 0*) \[Theta][0] == 1, \[Phi][0] == 1,(*Boundary Conditions as \[Eta]\[RightArrow]\[Infinity], approximated at \[Eta]= 50*) \[Theta][\[Eta]max] == 0, \[Phi][\[Eta]max] == 0}, { \[Theta], \[Phi]}, {\[Eta], 0, \[Eta]max}]; Join[{ ff} , sol\[Theta]\[Phi]]) & , {Function[\[Eta], -f0 Sqrt[1 + \[Eta]^2 ]], Function[\[Eta], Exp[-\[Eta] ]], Function[\[Eta], Exp[-\[Eta] ]]}, 3]; 

NestList converges well

Plot[{Through[solUN[[-2]][\[Eta] ]],Through[solUN[[-1]][\[Eta] ]]}, {\[Eta],0, \[Eta]max},PlotStyle -> {Automatic , {Dashed, Red}}] 

These results agree well with @bbgodfrey nice approach:

Plot[Through[solUN[[-1 ]][\[Eta] ]][[1]] - \[Eta], {\[Eta],0, \[Eta]max}, PlotRange -> All, AxesLabel -> {\[Eta], g[\[Eta]]}] 

Hope it helps!