**3. Analysis of Solution Stability**

By obeying the efforts of Merkin [61] and Merill et al. [62] from their outstanding discoveries of stability analysis scheme, the unsteady equations need to be deliberated in order to ultimately identify the reliable and stable solution since we notice the appearance of non-uniqueness solutions in the boundary value problem (7)–(9). Now, in accordance with the unsteady-state problem, a new similarity conversion is proposed

$$\begin{array}{l} \mu = \frac{dT}{1 - ct} \frac{\partial f}{\partial \eta}(\eta, \tau), \upsilon = -\left(\frac{dV}{1 - ct}\right)^{1/2} f(\eta, \tau), \theta(\eta, \tau) = \frac{T - T\_{\infty}}{T\_f - T\_{\infty}},\\ \eta = \sqrt{\frac{a}{\nu(1 - ct)}} y, \tau = \frac{a}{1 - ct} t. \end{array} \tag{14}$$

Employing the similarity variables of Equation (14) to Equations (7) and (8), we now obtain the following converted differential equations

$$\frac{\mu\_{\rm Imf}/\mu\_f}{\rho\_{\rm Imf}/\rho\_f}\frac{\partial^3 f}{\partial \eta^3} + \left(f + \frac{\varepsilon}{2}\eta\right)\frac{\partial^2 f}{\partial \eta^2} - \left(\frac{\partial f}{\partial \eta}\right)^2 - \varepsilon\frac{\partial f}{\partial \eta} - \left(1 + \varepsilon\tau\right)\frac{\partial^2 f}{\partial \eta\partial\tau} + \varepsilon + 1 = 0,\tag{15}$$

$$\frac{1}{\mathrm{Pr}} \frac{k\_{\mathrm{huf}}/k\_f}{\left(\rho \mathbb{C}\_p\right)\_{\mathrm{huf}}/\left(\rho \mathbb{C}\_p\right)\_f} \frac{\partial^2 \theta}{\partial \eta^2} + f \frac{\partial \theta}{\partial \eta} - 2\theta \frac{\partial f}{\partial \eta} - \frac{\varepsilon}{2} \eta \frac{\partial \theta}{\partial \eta} - \frac{\varepsilon}{2} 3\theta - \left(1 + \varepsilon \tau\right) \frac{\partial \theta}{\partial \tau} = 0,\tag{16}$$

with respect to 
$$f(0,\tau) = 0,\\ \frac{\partial f}{\partial \eta}(0,\tau) = \lambda + \gamma \frac{\partial^2 f}{\partial \eta^2}(0,\tau),\\ -\frac{k\_{\text{n}f}}{k\_f} \frac{\partial \mathcal{J}}{\partial \eta}(0,\tau) = \text{Bi}[1 - \theta(0,\tau)],\\ \text{as} \tag{17}$$

$$\frac{\partial f}{\partial \eta}(\eta,\tau) \to 1,\\ \theta(\eta,\tau) \to 0, \text{ as } \eta \to \infty.$$

In accordance with Weidman et al. [63], to test the stability of the steady flow *f*(η) = *f*0(η) and <sup>θ</sup>(η) = <sup>θ</sup>0(η) which fulfil the boundary value problem and boundary conditions (refer to (7)–(9)), we write

$$f(\eta, \tau) = f\mathfrak{o}(\eta) + e^{-\omega\tau}F(\eta), \quad \theta(\eta, \tau) = \theta\mathfrak{o}(\eta) + e^{-\omega\tau}G(\eta), \tag{18}$$

by which ω is the eigenvalue of unidentified variables, while functions *<sup>F</sup>*(η) and *<sup>G</sup>*(η) are relatively small to *f*0(η) and <sup>θ</sup>0(η). The eigenvalue problems (15) and (16) result in an infinite group of eigenvalues ω1 < ω2 < ω3...... that detect an early decay when ω1 is positive, while an early growth of disruptions

*Mathematics* **2020**, *8*, 1649

is observed when ω1 is negative, which exposes the unstable flow. Substituting (18) into (15)–(17), we develop

$$\frac{\mu\_{\rm Imf} / \mu\_f}{\rho\_{\rm Imf} / \rho\_f} \frac{\partial^3 F}{\partial \eta^3} + \left( f\_0 + \frac{\varepsilon}{2} \eta \right) \frac{\partial^2 F}{\partial \eta^2} + F \frac{\partial^2 f\_0}{\partial \eta^2} - 2 \frac{\partial f\_0}{\partial \eta} \frac{\partial F}{\partial \eta} + (\omega - \varepsilon) \frac{\partial F}{\partial \eta} = 0,\tag{19}$$

$$\frac{1}{\Pr} \frac{k\_{\rm mf} / k\_f}{\left(\rho \mathbb{C}\_p\right)\_{\rm mf} / \left(\rho \mathbb{C}\_p\right)\_f} \frac{\partial^2 G}{\partial \eta^2} + \left(f\_0 - \frac{\varepsilon}{2}\eta\right) \frac{\partial G}{\partial \eta} - 2\left(\theta\_0 \frac{\partial F}{\partial \eta} + G \frac{\partial f\_0}{\partial \eta}\right) + F \frac{\partial \theta\_0}{\partial \eta} + \left(\omega - \frac{3}{2}\varepsilon\right) G = 0,\tag{20}$$

and the boundary conditions are

$$\begin{split} F(0,\tau) = 0, & \frac{\partial F}{\partial \eta}(0,\tau) - \gamma \frac{\partial^2 F}{\partial \eta^2}(0,\tau) = 0, & -\frac{k\_{\rm mf}}{k\_f} \frac{\partial G}{\partial \eta}(0,\tau) - \text{Bi}G(0,\tau) = 0, \\ & \frac{\partial F}{\partial \eta}(\eta,\tau) \to 0, & G(\eta,\tau) \to 0, \text{ as } \eta \to \infty. \end{split} \tag{21}$$

The heat transfer stability and steady-state flow solutions *f*0(η) and <sup>θ</sup>0(η) was implemented via τ → 0, therefore *F* = *<sup>F</sup>*0(η) and *G* = *<sup>G</sup>*0(η) in (19)–(21). As a consequence, an early growth of Equation (18) is detected, and the subsequent generalized eigenvalue problem is recognized

$$\frac{\mu\_{\text{Im}f}/\mu\_f}{\rho\_{\text{Im}f}/\rho\_f} F\_0^{\prime\prime} + \left(f\_0 + \frac{\varepsilon}{2}\eta\right) F\_0^{\prime\prime} + F\_0 f\_0^{\prime\prime} - \left(2f\_0^{\prime} - \omega + \varepsilon\right) F\_0^{\prime} = 0,\tag{22}$$

$$\frac{1}{\mathrm{Pr}} \frac{k\_{\mathrm{lnf}}/k\_f}{\left(\rho \mathbb{C}\_p\right)\_{\mathrm{lnf}}/\left(\rho \mathbb{C}\_p\right)\_f} G\_0^{\prime\prime} + \left(f\mathrm{o} - \frac{\varepsilon}{2}\eta\right) G\_0^{\prime} + F\_0 \mathcal{O}\_0^{\prime} - 2\left(\mathcal{O}\_0 \mathrm{F}\_0^{\prime} + \mathrm{G}\_0 \mathrm{f}\_0^{\prime}\right) + \left(\omega - \frac{3}{2}\varepsilon\right) G\_0 = 0,\tag{23}$$

subject to

$$\begin{aligned} F\_0(0) &= 0, \newline F\_0'(0) - \gamma F\_0''(0) = 0, \newline -\frac{k\_{\rm{inf}}}{k\_f} G\_0'(0) - \text{Bi} G\_0(0) = 0, \\\ &F\_0'(\eta) \to 0, \newline G\_0(\eta) \to 0, \text{ as } \eta \to \infty. \end{aligned} \tag{24}$$

The range of possible eigenvalues can be calculated by resting a boundary condition [64] of the present problem. In this study, we choose to repose *F* 0(η) → 0, and the linear eigenvalue problems (22)–(24) are disclosed as *F* 0(0) = 1 for a fixed value of ω1. It is worth mentioning that the values of ω1 are proficient in measuring the stability of the corresponding solutions *f*0(η) and <sup>θ</sup>0(η).
