\[f'=y_2\\ f''=y_2'=y_3\\ f'''=y_3'=\frac{y_3y_5}{y_4-\theta_r}+\frac{y_4-\theta_r}{\theta_r}y_1y_3+\frac{y_4-\theta_r}{\theta_r}y_2^2-+\frac{y_4-\theta_r}{\theta_r}\left(\frac{2\beta}{(\eta+\alpha)^4}\right)y_4-\frac{y_4-\theta_r}{\theta_r}My_2\\ \theta'=y_5\\ \theta''=y_5'=-\frac{ay_5^2}{1+ay_4}-Pr\frac{y_1y_5}{1+ay_4}+\frac{2\beta\lambda(y_4+\epsilon)y_1}{(1+ay_4)(\eta+\alpha)^3}\]
Boundary conditions are:
\[y_1(0)=s,\ y_2(0)=1, \ y_4(0)=1, \ y_2(\infty)=0,\ y_4(\infty)=0\]
Set of equation (19) as well as boundary condition (20) are integrated numerically as an initial value problem to a given terminal point.