## 6 Numerical Scheme  
The system with high nonlinearity comprising Equations (12)-(15) with the support of Equation (16) is numerically solved by bvp4c MATLAB function. The following code transforms the given model into the 1 st order system of ODEs. \[ \begin{array}{l}g=y_{1},\\ g^{\prime}=y_{2},\\ g^{\prime\prime}=y_{3},\\ g^{\prime\prime\prime}=yy_{1},\\ yy_{1}=\frac{1}{2\xi}\Bigg{[}\begin{array}{l}\lambda y_{2}-2y_{3}+(1-\varphi )^{2.5}(1-\varphi+\varphi\frac{p_{CNT}}{\rho_{f}})F_{r}y_{2}^{2}\\ -(1-\varphi)^{2.5}\alpha\beta(y_{5}-y_{2})-(1-\varphi_{d})(1-\varphi)^{2.5}(1- \varphi+\varphi\frac{p_{CNT}}{\rho_{f}})y_{1}y_{3}\\ +(1-\varphi)^{2.5}\frac{M}{1+m^{2}y_{2}}\end{array}\Bigg{]},\\ \theta=y_{6},\\ \theta^{\prime}=y_{7},\\ \theta^{\prime\prime}=yy_{3},\end{array}  \tag{31}\]  
\[yy_{3}=\frac{1}{\frac{k_{sf}}{k_{sf}}\xi+\frac{4}{3N_{T}}(1+(\theta_{r}-1) \theta)^{3}}\left[ \begin{array}{l}\frac{k_{sf}}{k_{sf}}\xi y_{7}-0.5\text{Pr} (1-\varphi)^{2.5}(1-\varphi+\varphi\frac{p_{CNT}}{\rho_{f}})y_{1}y_{7}\\ -\frac{4}{N_{T}}(1+(\theta_{r}-1)\theta)^{2}\xi(\theta_{r}-1)\theta^{\prime 2 }\\ -\frac{4}{3N_{T}}(1+(\theta_{r}-1)\theta)^{3}0.5\theta-\frac{4E_{\text{Pr}}}{(1 -\varphi)^{2.5}}\xi y_{3}^{2}\\ -2\alpha\beta EcPr(y_{5}-y_{2})^{2}-2 EcPr(\frac{\lambda}{(1-\varphi)^{2.5}}+ \frac{M}{1+m^{2})y_{2}^{2}}\\ -0.5\alpha\beta_{T}\text{Pr}(1-\varphi+\varphi\frac{(p_{CNT})}{(C_{p})_{bf}})( y_{8}-y_{6})\end{array} \right] \tag{32}\]\[ \begin{split}&\theta_{p}=y_{8},\\ &\theta_{p}\prime=yy_{4},\\ & yy_{4}=\frac{\gamma\beta_{T}(1-\varphi+\varphi\frac{(C_{p})_{CNT }}{(C_{p})_{bf}})}{y_{4}}(y_{8}-y_{6}).\end{split}  \tag{34}\]  
With boundary conditions  
\[ \begin{split}& y_{1}(a)=\frac{a}{2}\varepsilon,y_{2}(a)=\frac{ \varepsilon}{2},y_{4}(a)=\frac{a}{2}\varepsilon,y_{5}(a)=\frac{\varepsilon}{ 2},\\ & y_{6}(a)=1,y_{2}(\infty)=0,y_{6}(\infty)=0,y_{8}(\infty)=0. \end{split}  \tag{35}\]  
Table 3 depicts the validation of the obtained results by comparing with already published articles in limiting case. This endorses the truthfulness of the presented mathematical model. Table 4 illustrates the numerically calculated values of the skin friction coefficient for numerous estimates of \(\alpha\), \(\phi_{d}\), \(F_{r}\), \(\beta\), \(M\), and \(m\). It is noticed that value of the drag force coefficient is higher for \(\alpha\), \(\phi_{d}\), \(F_{r}\), and \(M\), but it declines for estimates of \(\beta\) and \(m\) for both types of CNTs. Likewise, Table 5 portrays the Nusselt number for \(F_{r}\), \(\alpha\), \(\beta\), \(m\), \(N_{r}\), and \(\theta_{r}\). It is comprehended that heat transfer rate is higher in case of \(\alpha\) and \(F_{r}\), but converse behavior is seen for \(\beta\), \(m\), \(N_{r}\), and \(\theta_{r}\) for SWCNTs and MWCNTs. \begin{table}
\begin{tabular}{c c c c c c} \hline \hline \(a\) & \(\phi_{d}\) & \(F_{r}\) & \(\beta\) & \(M\) & \(m\) &
\begin{tabular}{c} **Skin Friction Coefficient** \\ **SWCNT** \\ \end{tabular} \\ \hline
0.001 & & & & & 0.00184647 & 0.00184157 \\
0.01 & & & & & 0.00583012 & 0.00581468 \\
0.2 & & & & & 0.02591430 & 0.02584880 \\  & 0.1 & & & & 0.00583012 & 0.00581468 \\  & 2.0 & & & & 0.00583276 & 0.00581722 \\  & 3.5 & & & & 0.00583491 & 0.00581927 \\  & 0.10 & & & & 0.00583012 & 0.00581468 \\  & 0.25 & & & & 0.00644651 & 0.00640784 \\  & 0.4 & & & & 0.00706302 & 0.00700110 \\  & 1.0 & & & 0.00583012 & 0.00581468 \\  & 2.0 & & & 0.00582586 & 0.00581043 \\  & 3.0 & & & 0.00582165 & 0.00580623 \\  & 0.2 & & & 0.00583012 & 0.00581468 \\  & 0.3 & & & 0.00711701 & 0.00710157 \\  & 0.4 & & & 0.00840431 & 0.008388888 \\  & 1.0 & & 0.00583012 & 0.00581468 \\  & 1.4 & & 0.00499559 & 0.00498015 \\  & 1.8 & & 0.00447082 & 0.00445537 \\ \hline \hline \end{tabular}
\end{table} 
Table 4: Skin friction coefficient against different parameters.