**5. Stability Analysis**

Merkin [39,40] established an improved version of the stability analysis, which is prominent among researchers for the examination of the stability of numerical solutions. Because we observed dual solutions in the present work, we assess the solution's stability to determine the flow behavior. To initiate the linear stability analysis, the model equations in Equations (3)–(5) need to be considered in the unsteady form as follows:

$$\frac{\partial \overline{u}}{\partial \overline{x}} + \frac{\partial \overline{v}}{\partial \overline{y}} = 0,\tag{28}$$

$$\frac{\partial\overline{u}}{\partial\overline{t}} + \overline{u}\frac{\partial\overline{u}}{\partial\overline{x}} + \overline{v}\frac{\partial\overline{u}}{\partial\overline{y}} = \frac{\mu\_{nf}}{\rho\_{nf}}\frac{\partial^2\overline{u}}{\partial\overline{y}^2} + \frac{\mu\_{nf}}{\rho\_{nf}}\frac{m\mathbb{E}^2}{2}\left(\frac{\partial\overline{u}}{\partial\overline{y}}\right)^2\frac{\partial^2\overline{u}}{\partial\overline{y}^2} + \overline{u}\_\varepsilon\frac{d\overline{u}\_\varepsilon}{d\overline{x}} - \frac{\sigma B\_0^2(\overline{u}\_\varepsilon - \overline{u})}{\rho\_{nf}},\tag{29}$$

$$\frac{\partial T}{\partial \overline{t}} + \overline{u}\frac{\partial T}{\partial \overline{\mathbf{x}}} + \overline{v}\frac{\partial T}{\partial \overline{y}} = \frac{k\_{nf}}{\left(\rho \mathbb{C}\_p\right)\_{nf}}\frac{\partial^2 T}{\partial \overline{y}^2} + \frac{16\,\sigma\_1 T\_{\text{os}}^3}{3\left(\rho \mathbb{C}\_p\right)\_{nf}k\_1}\frac{\partial^2 T}{\partial \overline{y}^2} \tag{30}$$

with the boundary conditions of Equation (6). Then, we introduce the dimensional time variable, *t* with a new similarity variable τ = τ0/*x*4/5 through scaling group analysis. The new similarity transformation is given as:

$$\begin{aligned} \eta &= \frac{y}{x^5}, \quad \psi = x^{\frac{3}{5}} f(\eta, \tau), \quad \sigma = \sigma\_0 \ge^{-\frac{4}{5}}, \quad \theta = \theta\_0(\eta, \tau) \ge^{\frac{2}{5}},\\ \mu\_t &= (u\_t)\_0 x^{\frac{1}{5}}, \quad u\_1 = (u\_1)\_0 \ge^{-\frac{1}{5}}, \quad m = m\_0 \ge^{\frac{2}{5}}, \quad \tau = \frac{\tau\_0}{x^{\frac{4}{5}}}. \end{aligned} \tag{31}$$

Employing (31) in the dimensionless form of Equations (28)–(30) and (6) gives the following system of equations:

$$\frac{A\_1}{A\_2} \frac{\partial^3 f}{\partial \eta^3} \left[ 1 + \frac{m\_0 D \varepsilon}{2} \left( \frac{\partial^2 f}{\partial \eta^2} \right)^2 \right] - \frac{1}{5} \left( \frac{\partial f}{\partial \eta} \right)^2 + \frac{3}{5} f(\eta, \tau) \frac{\partial^2 f}{\partial \eta^2} + \frac{1}{5} - \frac{M}{A\_2} \left( 1 - \frac{\partial f}{\partial \eta} \right) - \frac{\partial^2 f}{\partial \eta \partial \tau} = 0,\tag{32}$$

$$\left(A\_4 + \frac{4}{3}Rd\right)\frac{\partial^2 \theta}{\partial \eta^2} - \frac{2}{5}A\_3 \text{Pr}\frac{\partial f}{\partial \eta}\,\theta(\eta, \tau) + \frac{3}{5}A\_3 \text{Pr}\,f(\eta, \tau)\,\frac{\partial \theta}{\partial \eta} - \frac{\partial \theta}{\partial \tau} = 0,\tag{33}$$

with the boundary conditions:

$$f(0,\tau) = \frac{5}{3} f\_{\text{uv}}, \quad \frac{\partial f}{\partial \eta}(0,\tau) = \varepsilon, \quad \theta(0,\tau) = 1, \quad \frac{\partial f}{\partial \eta}(\infty,\tau) = 1, \quad \theta(\infty,\tau) = 0. \tag{34}$$

It is assumed that the solutions of (32)–(34) are expressed by the formulas of Equation (35):

$$f(\eta, \tau) = f\_0(\eta) + \varepsilon^{-\gamma \tau} F(\eta, \tau), \quad \theta(\eta, \tau) = \theta\_0(\eta) + \varepsilon^{-\gamma \tau} G(\eta, \tau), \tag{35}$$

where *f*(η) = *f*0(η) and <sup>θ</sup>(η) = <sup>θ</sup>0(η) are the solutions found in the previous section, in which the disturbance is superimposed to determine their stability. Here, the unknown eigenvalue parameter is denoted by γ, and *<sup>F</sup>*(η, τ) and *<sup>G</sup>*(η, τ) are relatively small compared to the steady state solutions (*f*0(η) and <sup>θ</sup>0(η)). The substitution of Equation (35) into Equations (32)–(34) gives the following system:

$$\frac{\partial^3 F}{\partial \eta^3} \Big( 1 + \frac{A\_1}{A\_2} \frac{m\_0 \mathcal{D} e}{2} \left. f\_0'' \right|^2 \Big) + \frac{\partial^2 F}{\partial \eta^2} \Big( \frac{3}{5} f\_0 - \frac{2}{5} f\_0' \Big) + \frac{3}{5} f\_0'' F + \left( \frac{M}{A\_2} + \gamma \right) \frac{\partial F}{\partial \eta} - \frac{\partial^2 F}{\partial \eta \partial \tau} = 0,\tag{36}$$

$$\left(A\_4 + \frac{4}{3}Rd\right)\frac{\partial^2 G}{\partial \eta^2} + \frac{3}{5}A\_3 \text{Pr} f\_0 \frac{\partial G}{\partial \eta} + \frac{3}{5}A\_3 \text{Pr} F\theta\_0 - \frac{2}{5}A\_3 \text{Pr} f\_0' G - \frac{2}{5}A\_3 \text{Pr} \theta\_0 \frac{\partial F}{\partial \eta} + \gamma G - \frac{\partial G}{\partial \tau} = 0,\tag{37}$$

subject to the boundary conditions:

$$F(0, \tau) = 0, \quad \frac{\partial F}{\partial \eta}(0, \tau) = 0, \quad G(0, \tau) = 0, \quad \frac{\partial F}{\partial \eta}(\infty, \tau) = 0, \quad G(\infty, \tau) = 0. \tag{38}$$

Referring to Merkin [39,40], τ → 0 is fixed to examine the stability of the steady state boundary layer flow. Thus, *F* = *<sup>F</sup>*0(η) and *G* = *<sup>G</sup>*0(η) in Equations (37)–(39), yielding the following linearized eigenvalue problem:

$$F\_{0}^{\prime\prime}\left(1+\frac{A\_{1}}{A\_{2}}\frac{m\_{0}D\varepsilon}{2}\left(f\_{0}^{\prime\prime}\right)^{2}\right) + F\_{0}^{\prime}\left(\frac{3}{5}f\_{0}-\frac{2}{5}f\_{0}^{\prime}\right) + \frac{3}{5}f\_{0}^{\prime\prime}F\_{0} + \left(\frac{M}{A\_{2}}+\gamma\right)F\_{0}^{\prime} = 0,\tag{39}$$

$$\left(A\_4 + \frac{4}{3}Rd\right)\mathbf{G}\_0' + \frac{3}{5}A\_3 \mathbf{Pr} f\_0 \mathbf{G}\_0' + \frac{3}{5}A\_3 \mathbf{Pr} F\_0 \theta\_0 - \frac{2}{5}A\_3 \mathbf{Pr} f\_0' \mathbf{G}\_0 - \frac{2}{5}A\_3 \mathbf{Pr} \theta\_0 \mathbf{F}\_0' + \gamma \mathbf{G}\_0 = 0,\tag{40}$$

with the boundary conditions:

$$\begin{array}{llll} F\_0(\eta) = 0, & F\_0'(\eta) = 0, & G\_0(\eta) = 0, & \text{at} & \eta = 0, \\ F\_0'(\eta) = 0, & G\_0(\eta) = 0 & \text{as} & \eta \to \infty. \end{array} \tag{41}$$

It is necessary to replace one of the outer boundary conditions with a normalizing boundary condition to obtain the eigenvalues. Therefore, the boundary condition *<sup>F</sup>*0(∞) = 0 is substituted with *F*0 (0) = 1. The system of equations in Equations (38)–(40) with the new boundary condition is solved by the MATLAB boundary value problem solver (bvp4c) to obtain the lowest eigenvalues as the governing parameter varies. These lowest eigenvalues are classified according to their sign. If the lowest eigenvalue falls in the positive range of values, then the respective numerical solution is accepted as a stable solution. Meanwhile, the negative lowest eigenvalue suggests the numerical solution is unstable. Further explanation about the stable and unstable solutions is provided in the next section.
