## 3 Stability Analysis  
The results display the presence of dual solutions for the tested parameter. Due to the importance of identifying the physically-realizable solution, a stability analysis is performed. Following Merkin [30], a new dimensionless variable, \(\tau\), is introduced, which is related to an initial value problem. The new dimensionless variables can be represented in Equation (18):  
\[ \begin{array}{l}\eta=\Big{(}\frac{a}{2 lv_{f}(1-ct)}\Big{)}^{1/2}e^{x/2l}y,\; \psi=\Big{(}\frac{2 lv_{p,\mu}}{1-ct}\Big{)}^{1/2}e^{x/2l}f(\eta,\tau),\;\theta (\eta,\tau)=\frac{T-T_{\infty}}{T_{w}-T_{\infty}},\\ \tau=\frac{at}{2l(1-ct)}e^{x/l}.\end{array}  \tag{18}\]  
Equation (18) is then substituted into (2) and (3) to obtain  
\[ \begin{array}{l}\frac{1}{\left(1-\varphi\right)^{2.5}\Big{(}1-\varphi+ \varphi\Big{(}\frac{c_{p}}{F_{f}}\Big{)}\Big{)}}\frac{\partial^{3}f}{\partial y ^{3}}+f\frac{\partial^{2}f}{\partial\eta^{2}}-2\Big{(}\frac{\partial f}{ \partial\eta}\Big{)}^{2}+2-A\Big{(}2\frac{\partial f}{\partial\eta}+\eta \frac{\partial^{2}f}{\partial\eta^{2}}-2\Big{)}\\ -\Big{(}2\tau\frac{\partial f}{\partial\eta}+2A\tau+1\Big{)}\frac{\partial^{2} f}{\partial\eta\sigma\tau}=0,\end{array}  \tag{19}\]  
\[\frac{1}{\Pr}\frac{k_{nf}}{k_{f}\Big{(}1-\varphi+\varphi\Big{(}\frac{(\rho c _{p})_{x}}{\left(\rho c_{p}\right)_{f}}\Big{)}\Big{)}}\frac{\partial^{2} \theta}{\partial\eta^{2}}+f\frac{\partial\theta}{\partial\eta}-\theta\frac{ \partial f}{\partial\eta}-A\bigg{(}2\theta+\eta\frac{\partial\theta}{\partial \eta}\bigg{)}-(1+2A\tau+2\tau)\frac{\partial\theta}{\partial\tau}=0, \tag{20}\]  
along with the boundary conditions:  
\[ \begin{array}{l}f(0,\tau)=s,\;\frac{\partial f}{\partial\eta}(0,\tau)=\varepsilon +\sigma\frac{\partial^{2}f}{\partial\eta^{2}}(0,\tau),\;\theta(0,\tau)=1,\\ \frac{\partial f}{\partial\eta}(\eta,\tau)\to 1,\;\theta(\eta,\tau)\to 0,\;\text{as}\;\eta \rightarrow\infty.\end{array}  \tag{21}\]The stability of solutions \(f(\eta)=f_{0}(\eta)\) and \(\theta(\eta)=\theta_{0}(\eta)\) that satisfy the boundary value problem Equations (11)-(13) is tested by introducing (see [PERSON] [30] and [PERSON] et al. [31]):  
\[f(\eta,\tau)=f_{0}(\eta)+e^{-\,\gamma\,\tau}F(\eta,\tau)\text{ and }\theta(\eta, \tau)=\theta_{0}(\eta)+e^{-\,\gamma\,\tau}G(\eta,\tau), \tag{22}\]  
where \(\gamma\) is an unknown eigenvalue parameter, \(F(\eta,\tau)\) and \(G(\eta,\tau)\) are small corresponding to \(f_{0}(\eta)\) and \(\theta_{0}(\eta)\), respectively. The linearized problem can be obtained by substituting Equation (22) into (19) and (20), as below:  
\[ \begin{split}&\frac{1}{(1-\varphi)^{2.5}\left(1-\varphi+\varphi \left(\frac{\rho_{\tau}}{\rho_{f}}\right)\right)}\frac{\partial^{2}F}{ \partial\eta^{3}}+f_{0}\frac{\partial^{2}F}{\partial\eta^{2}}+Ff_{0}{}^{ \prime\prime}-4f_{0}{}^{\prime}\frac{\partial F}{\partial\eta}-A\eta\frac{ \partial^{2}F}{\partial\eta^{2}}\\ &-2A\frac{\partial F}{\partial\eta}-2A\frac{\partial^{2}F}{ \partial\eta\partial\tau}+2A\gamma\frac{\partial F}{\partial\eta}\frac{ \partial^{2}F}{\partial\eta\tau}+\gamma\frac{\partial F}{\partial\eta}-2\tau \frac{\partial^{2}F}{\partial\eta\partial\tau}+2\tau\gamma\frac{\partial F}{ \partial\eta}=0,\\ &\frac{1}{\Pr}\frac{k_{nf}}{k_{f}\left(1-\varphi+\varphi\left( \frac{\rho_{\tau}}{\rho_{f}}\right)\right)}\frac{\partial^{2}G}{\partial\eta ^{2}}-A\eta\frac{\partial G}{\partial\eta}-2A\tau\frac{\partial G}{\partial \tau}-2A\tau\gamma G-\frac{\partial G}{\partial\tau}+\gamma G\\ &-2 AG-2\tau f_{0}{}^{\prime}\frac{\partial G}{\partial\tau}+\gamma 2 \tau f_{0}{}^{\prime}G-\theta_{0}\frac{\partial F}{\partial\eta}-f_{0}{}^{ \prime}G+\theta_{0}{}^{\prime}F+f_{0}\frac{\partial G}{\partial\eta}=0,\end{split}  \tag{23}\]  
subject to boundary conditions  
\[ \begin{split}& F(0,\tau)=0,\ \frac{\partial F(0,\tau)}{ \partial\eta}+\sigma\frac{\partial^{2}F(0,\tau)}{\partial\eta^{2}}=0,\ G(0, \tau)=0,\\ &\frac{\partial F(\eta,\tau)}{\partial\eta}\to 0,\ G(\eta,\tau) \to 0\ \ \ \text{as }\eta\rightarrow\infty.\end{split}