**4. Implicit Finite Difference Method**

The implicit finite difference method (I-FDM) [40] was utilized to obtain the numerical solution for the equation set of models. I-FDM has an advantage as it has fast convergence. It is 2nd order convergen<sup>t</sup> and inherently stable. I-FDM satisfies the Von Neumann stability test, which has the criterion of a real numerical solution for PDEs with the help of the stability and consistency of a numerical solution. I-FDM is employed to obtain the solution of Equations (11) and (12) using boundary conditions (13). This is a suitable method to obtain the approximated solution of boundary layer problems. I-FDM is widely applicable in the flow problems of the laminar boundary layer, and the obtained results are more effective than others.

To apply the implicit finite difference method [41], Equations (11) and (12) were written in the form of 1st order differential equations utilizing newly employed variables. Reduced equations are as follows [32]:

$$L\_1 = f',\tag{16}$$

$$L\_2 = L\_{1\prime}^{\prime} \tag{17}$$

$$z\_3 = \theta',\tag{18}$$

$$\tau^\* L\_2' \left(1 - \varsigma^\* L\_2 \right)^2 + \phi\_{\bar{\mathbf{Y}}\_2} \left[f L\_2 - L\_1^2 \right] - \frac{1}{\phi\_{\bar{\mathbf{Y}}\_1}} F\_\pi L\_1 = 0,\tag{19}$$

$$L\_3' \left(1 + \nu^\* \theta + \frac{1}{\Phi \dot{\chi}\_4} P\_r N\_\pi \right) + \nu^\* L\_3^2 + P\_r \frac{\phi\_{\mathcal{Y}\_3}}{\Phi \dot{\chi}\_4} [f L\_3 - L\_1 \theta] = 0. \tag{20}$$

With the presence of newly employed variables, boundary conditions eventually changed to [31]

$$f(0) = S,\\ L\_1(0) = 1 + \Lambda\_\pi \, L\_2(0),\\ L\_3(0) = -B\_\pi(1 - \theta(0)),\\ L\_1(\infty) \to 0,\\ \theta(\infty) \to 0. \tag{21}$$

The different formulae were calculated using central differencing, and average functions were replaced. Thus, the 1st ODEs (16) and (20) order decreases to the next series of nonlinear algebraic formulae.

$$\frac{(L\_1)\_j + (L\_1)\_{j-1}}{2} = \frac{f\_j - f\_{j-1}}{h},\tag{22}$$

$$\frac{(L\_2)\_j + (L\_2)\_{j-1}}{2} = \frac{(L\_1)\_j - (L\_1)\_{j-1}}{h},\tag{23}$$

$$\frac{(L\_3)\_j + (L\_3)\_{j-1}}{2} = \frac{\theta\_j - \theta\_{j-1}}{h},\tag{24}$$

$$\begin{split} &\pi^\* \left( \frac{(l\_{\Sigma})\_{j} - (l\_{\Sigma})\_{j-1}}{h} \right) \left( 1 - \xi^\* \left( \frac{(l\_{\Sigma})\_{j} + (l\_{\Sigma})\_{j-1}}{2} \right)^2 \right) \\ &+ \left[ \Phi \chi\_{2} \left( \left( \frac{f\_{j} + f\_{j-1}}{2} \right) \left( \frac{(l\_{\Sigma})\_{j} + (l\_{\Sigma})\_{j-1}}{2} \right) - \left( \frac{(l\_{\Sigma})\_{j} + (l\_{\Sigma})\_{j-1}}{2} \right)^2 \right) - F\_{\pi} \frac{1}{\Phi \chi\_{1}} \left( \frac{(l\_{\Sigma})\_{j} + (l\_{\Lambda})\_{j-1}}{2} \right) \right], \end{split} \tag{25}$$

$$\begin{split} & \left( \frac{(L\_3)\_j - (L\_3)\_{j-1}}{h} \right) \left( 1 + \upsilon^\* \left( \frac{\theta\_j + \theta\_{j-1}}{2} \right) + \frac{1}{\Phi\_{\overline{\chi}\_4}} P\_r N\_{\overline{\pi}} \right) + \upsilon^\* \left( \frac{(L\_3)\_j + (L\_3)\_{j-1}}{2} \right)^2 \\ & + P\_r \frac{\Phi\_{\overline{\chi}\_3}}{\Phi\_{\overline{\chi}\_4}} \left[ \left( \frac{f\_j + f\_{j-1}}{2} \right) \left( \frac{(L\_3)\_j + (L\_3)\_{j-1}}{2} \right) - \left( \frac{(L\_1)\_j + (L\_1)\_{j-1}}{2} \right) \left( \frac{\theta\_j + \theta\_{j-1}}{2} \right) \right] = 0. \end{split} \tag{26}$$

To linearize the resulting equations, Newton's technique was used. As an example, consider iteration (*i* + 1)th

$$\mathbf{O}(\mathbf{j}\_j^{(i+1)}=\mathbf{O}\_j^{(i)}+\mathbf{O}(\mathbf{j}\_j^{(i)}.\tag{27}$$

under the substitutition of the linear tridiagonal equational system into Equations (22)–(26), disregarding the elevated O ¨ *i j* components.

$$\bullet f\_{\hat{\jmath}} - \bullet f\_{\hat{\jmath}-1} - \frac{1}{2}h(\bullet (L\_1)\_{\hat{\jmath}} + \bullet (L\_1)\_{\hat{\jmath}-1}) = (d\_1)\_{\hat{\jmath}-\frac{1}{2}\prime} \tag{28}$$

$$\dot{\mathcal{O}}(L\_1)\_{\dot{j}} - \ddot{\mathcal{O}}(L\_1)\_{\dot{j}-1} - \frac{1}{2}h(\ddot{\mathcal{O}}(L\_2)\_{\dot{j}} + \ddot{\mathcal{O}}(L\_2)\_{\dot{j}-1}) = (d\mathbf{2})\_{\dot{j}-\frac{1}{2}\prime} \tag{29}$$

$$\ddot{\mathcal{O}}\theta\_{\dot{\jmath}} - \ddot{\mathcal{O}}\theta\_{\dot{\jmath}-1} - \frac{1}{2}h(\ddot{\mathcal{O}}(L\_{\eth})\_{\dot{\jmath}} + \ddot{\mathcal{O}}(L\_{\eth})\_{\dot{\jmath}-1}) = (d\mathfrak{d})\_{\dot{\jmath}-\frac{1}{2}\prime} \tag{30}$$

$$\begin{aligned} (a\_1)\_j \bullet f\_j + (a\_2)\_j \bullet f\_{j-1} + (a\_3)\_j \bullet \mathcal{U}\_{1j} + (a\_4)\_j \bullet \mathcal{U}\_{j-1} + (a\_5)\_j \bullet \mathcal{U}\_{2j} + (a\_6)\_j \bullet \mathcal{U}\_{2j-1} \\ + (a\_7)\_j \bullet \mathcal{O}\_{\underline{j}} + (a\_8)\_{\underline{j}} \bullet \mathcal{O}\_{\underline{j}-1} + (a\_9)\_{\underline{j}} \bullet (\mathcal{U}\_{3j})\_{\underline{j}} + (a\_{10})\_{\underline{j}} \bullet (\mathcal{U}\_{3})\_{\underline{j}-1} = (d\_4)\_{j - \frac{1}{2}} \end{aligned} \tag{31}$$

$$\begin{split} \bullet \ (b\_1)\_j \ddot{\mathsf{O}} f\_j + (b\_2)\_j \ddot{\mathsf{O}} f\_{j-1} + (b\_3)\_j \ddot{\mathsf{O}} L\_{1j} + (b\_4)\_j \ddot{\mathsf{O}} L\_{1j-1} + (b\_5)\_j \ddot{\mathsf{O}} L\_{2j} + (b\_6)\_j \ddot{\mathsf{O}} L\_{2j-1} \\ \bullet + (b\gamma)\_j \ddot{\mathsf{O}} \theta\_j + (b\varsigma)\_j \ddot{\mathsf{O}} \theta\_{j-1} + (b\varsigma)\_j \ddot{\mathsf{O}} (L\_3)\_j + (b\_{10})\_j \ddot{\mathsf{O}} (L\_3)\_{j-1} = (d\mathfrak{s})\_{j-\frac{1}{2}}. \end{split} \tag{32}$$

where

$$(d\_1)\_{j-\frac{1}{2}} = -f\_j + f\_{j-1} + \frac{h}{2}(L\_1)\_j + ((L\_1)\_{j-1})\_j \tag{33}$$

$$(d\_2)\_{j - \frac{1}{2}} = -(L\_1)\_j + (L\_1)\_{j-1} + \frac{h}{2}((L\_2)\_j + (L\_2)\_{j-1}),\tag{34}$$

$$(d\_3)\_{j-\frac{1}{2}} = -\theta\_{\dot{\jmath}} + \theta\_{\dot{\jmath}-1} + \frac{h}{2}((L\_3)\_{\dot{\jmath}} + (L\_3)\_{\dot{\jmath}-1}),\tag{35}$$

(*d*4)*j*− 12 = −*<sup>h</sup>*⎡⎣*τ*<sup>∗</sup> (*<sup>L</sup>*2)*j* − (*<sup>L</sup>*2)*j*−<sup>1</sup> *h* ⎛⎝1 − *ς*∗ (*<sup>L</sup>*2)*j* + (*<sup>L</sup>*2)*j*−<sup>1</sup> 2 2⎞⎠⎤⎦ − *<sup>h</sup>*⎡⎣⎡⎣*<sup>φ</sup>b*⎛⎝ *fj* + *fj*−<sup>1</sup> 2 (*<sup>L</sup>*2)*j* + (*<sup>L</sup>*2)*j*−<sup>1</sup> 2 − (*<sup>L</sup>*1)*j* + (*<sup>L</sup>*1)*j*−<sup>1</sup> 2 2⎞⎠ − *Fπ* 1*φa* (*<sup>L</sup>*1)*j* + (*<sup>L</sup>*1)*j*−<sup>1</sup> 2 ⎤⎦⎤⎦, (36)

(*d*5)*j*− 12 = −*h* (*<sup>L</sup>*3)*j*<sup>−</sup>(*<sup>L</sup>*3)*j*−<sup>1</sup> *h* 1 + *υ*<sup>∗</sup>*<sup>θ</sup>j*+*θj*−<sup>1</sup> 2 + 1*φ*Υ¨ 4 *PrNπ* − *h υ*<sup>∗</sup>(*<sup>L</sup>*3)*j*+(*<sup>L</sup>*3)*j*−<sup>1</sup> 2 2 −*hPr φ*Υ¨ 3 *φ*Υ¨ 4 (*fj*+*fj*−<sup>1</sup>)((*<sup>L</sup>*3)*j*+(*<sup>L</sup>*3)*j*−<sup>1</sup>) 4 + *hPr φ*Υ¨ 3 *φ*Υ¨ 4 (*<sup>θ</sup>j*+*θj*−<sup>1</sup>)((*<sup>L</sup>*1)*j*+(*<sup>L</sup>*1)*j*−<sup>1</sup>) 4 (37)

The boundary conditions become

$$
\ddot{\mathcal{O}}f\_0 = 0, \ddot{\mathcal{O}}(z\_1)\_0 = 0, \ddot{\mathcal{O}}(z\_3)\_0 = 0, \ddot{\mathcal{O}}(z\_1)\_f = 0, \ddot{\mathcal{O}}\theta\_f = 0. \tag{38}
$$

The following are the formulae (33)–(37) that produce the bulk tridiagonal array,

$$R\ddot{\mathbf{O}} = p,\tag{39}$$

where

$$\mathbf{H} = \begin{bmatrix} \begin{array}{c} \mathcal{O}\_{1} & \mathcal{O}\_{1} \\ \mathcal{O}\_{2} & \mathcal{O}\_{2} & \mathcal{E}\_{2} \\ & \ddots & \ddots & \ddots \\ & & & \ddots & \ddots & \ddots \\ & & & & \mathcal{O}\_{j-1} & \mathcal{O}\_{j-1} \\ & & & & & \mathcal{O}\_{j-1} & \mathcal{O}\_{j-1} \end{array} \end{bmatrix}, \mathbf{O} = \begin{bmatrix} \mathcal{O}\_{1} \\ \dot{\mathcal{O}}\_{2} \\ \vdots \\ \dot{\mathcal{O}}\_{j-1} \\ \dot{\mathcal{O}}\_{j} \end{bmatrix}, p = \begin{bmatrix} (d\_{1})\_{j-\frac{1}{2}} \\ (d\_{2})\_{j-\frac{1}{2}} \\ \vdots \\ (d\_{f-1})\_{j-\frac{1}{2}} \\ (d\_{f})\_{j-\frac{1}{2}} \end{bmatrix}. \tag{40}$$

This matrix, *H*, resembles the generalized size of *J* × *J*, whereas the O and ¨ *p* indicate the column vectors order of *J* × 1. Afterward, a unique LU factorization approach was employed to ge<sup>t</sup> the solution for O. ¨
