## 3 Stability Analysis  
The present numerical results depict the existence of dual solutions in a certain range of \(\lambda\), \(N\), \(K\), and \(M\). Therefore, the stability of dual solutions has to be determined. In fact, prior studies have established that the stability of dual solutions can be determined by performing stability analysis, as suggested by [PERSON] [4]. Thus, for the purpose of conducting the stability analysis, an unsteady-state problem is considered. Equation (1) holds, while Equations (2)-(4) are substituted as follows:  
\[\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u }{\partial y}+w\frac{\partial u}{\partial z}=v_{nf}\frac{\partial^{2}u}{ \partial z^{2}}-\frac{\sigma_{nf}B^{2}u}{\rho_{nf}}+\frac{g(\rho\beta)_{nf}}{ \rho_{nf}}(T-T_{\infty}), \tag{18}\]  
\[\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v }{\partial y}+w\frac{\partial v}{\partial z}=v_{nf}\frac{\partial^{2}v}{ \partial z^{2}}-\frac{\sigma_{nf}B^{2}v}{\rho_{nf}}, \tag{19}\]\[\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{ \partial y}+w\frac{\partial T}{\partial z}=\alpha_{nf}\frac{\partial^{2}T}{ \partial z^{2}}, \tag{20}\]  
where \(t\) denotes the time. Additionally, the new dimensionless functions of \(u\), \(v\), \(w\), and \(\theta\), are initiated in relation to the similarity variable, \(\eta\) and the second similarity variable, \(\tau\) for the unsteady-state problem (18)-(20) as follows:  
\[ \begin{array}{l}u=a(x+y)^{n}\frac{\partial f}{\partial\eta}(\eta,\tau),\ \ \ v=a(x+y)^{n}\frac{\partial u}{\partial\eta}(\eta,\tau),\\ w=-\sqrt{uv_{f}}(x+y)^{(n-1)/2}\Big{[}\frac{n+1}{2}(f+h)+\frac{n-1}{2}\eta \Big{(}\frac{\partial f}{\partial\eta}+\frac{\partial u}{\partial\eta}\Big{)} +(n-1)\tau\Big{(}\frac{\partial f}{\partial\tau}+\frac{\partial u}{\partial \tau}\Big{)}\Big{]},\\ \theta(\eta,\tau)=\frac{T-T_{\infty}}{T_{\infty}-T_{\infty}},\ \ \eta=\sqrt{\frac{a}{v_{f}}}(x+y)^{(n-1)/2}z,\ \ \ \tau=a(x+y)^{n-1}t.\end{array}  \tag{21}\]  
Thus, the unsteady-state problem (18)-(20) can be written as:  
\[ \begin{array}{l}\frac{1}{(1-\phi)^{25}}\frac{\partial^{3}f}{\partial\eta^{3 }}+\Big{(}1-\phi+\phi\frac{\rho_{s}}{\rho_{f}}\Big{)}\Big{[}\frac{n+1}{2}(f+h )\frac{\partial^{2}f}{\partial\eta^{2}}-n\Big{(}\frac{\partial f}{\partial \eta}+\frac{\partial u}{\partial\eta}\Big{)}\frac{\partial f}{\partial\eta}\\ -(n-1)\tau\Big{(}\frac{\partial f}{\partial\eta}+\frac{\partial u}{\partial \eta}\Big{)}\frac{\partial^{2}f}{\partial\eta\sigma T}+(n-1)\tau\Big{(}\frac{ \partial f}{\partial\tau}+\frac{\partial u}{\partial\tau}\Big{)}\frac{\partial ^{2}f}{\partial\eta^{2}}-\frac{\partial^{2}f}{\partial\eta\sigma T}\Big{]}\\ -\Big{[}1+\frac{3\left(\sigma_{s}/\sigma_{f}-1\right)\phi}{\sigma_{s}/\sigma_ {f}+2-\left(\sigma_{s}/\sigma_{f}-1\right)\phi}\Big{]}M\frac{\partial f}{ \partial\eta}+\Big{(}1-\phi+\phi\frac{\rho(\rho\beta_{s})}{\rho_{f}}\Big{)} \lambda\theta=0,\end{array}  \tag{22}\]  
\[ \begin{array}{l}\frac{1}{(1-\phi)^{25}}\frac{\partial^{3}h}{\partial\eta^{3 }}+\Big{(}1-\phi+\phi\frac{\rho_{s}}{\rho_{f}}\Big{)}\Big{[}\frac{n+1}{2}(f+h )\frac{\partial^{2}h}{\partial\eta^{2}}-n\Big{(}\frac{\partial f}{\partial \eta}+\frac{\partial u}{\partial\eta}\Big{)}\frac{\partial u}{\partial\eta} \\ -(n-1)\tau\Big{(}\frac{\partial f}{\partial\eta}+\frac{\partial u}{\partial \eta}\Big{)}\frac{\partial^{2}f}{\partial\eta\sigma T}+(n-1)\tau\Big{(}\frac{ \partial f}{\partial\tau}+\frac{\partial u}{\partial\tau}\Big{)}\frac{ \partial^{2}h}{\partial\eta^{2}}-\frac{\partial^{2}h}{\partial\eta\sigma T} \Big{]}\\ -\Big{[}1+\frac{3\left(\sigma_{s}/\sigma_{f}-1\right)\phi}{\sigma_{s}/\sigma_ {f}+2-\left(\sigma_{s}/\sigma_{f}-1\right)\phi}\Big{]}M\frac{\partial h}{ \partial\eta}=0,\end{array}  \tag{23}\]  
\[