**4. Methodology**

The momentum and energy coupled non-linear ODEs (Equations (9) and (10)) and (Equations (11) and (12)), along with the boundary conditions (BCs) in Equation (13) and Equation (14), are solved numerically via bvp4c in MATLAB, which is based on a threestage Lobatto technique for the various comprising parameters and gamma nanofluids. The three-stage Lobatto technique is a collocation technique with fourth-order accuracy. The form of the ODEs (ordinary differential equations), along with the BCs, is altered into the group of first order IVP (intial value problem) by exercising the new variables. This process is carried forward by introducing the following variables:

$$\mathbf{F} = \mathbf{Z}\_1, \mathbf{F}' = \mathbf{Z}\_2, \mathbf{F}'' = \mathbf{Z}\_3, \boldsymbol{\theta} = \mathbf{Z}\_{4\prime}, \boldsymbol{\theta}' = \mathbf{Z}\_{5\prime} \tag{28}$$

Utilizing Equation (28) for the aforementioned ODEs and the boundary conditions, we get a system of ODEs for the model of (γAl2O3−H2O) and (γAl2O3−C2H6O2) nanofluids, respectively given as:

$$\begin{pmatrix} d \\ d \\ d\eta \\ Z\_{2} \\ Z\_{5} \end{pmatrix} = \begin{pmatrix} Z\_{2} \\ Z\_{3} \\ \hline \\ Z\_{5} \\ Z\_{6} \\ \hline \\ \begin{pmatrix} -4R\_{d}\mathbf{K}\_{6}(1+(\theta\_{w}-1)Z\_{4})^{2} - \mathbf{K}\_{5}\mathbf{Z}\_{2} + \mathbf{K}\_{4}\mathbf{Z}\_{4} \\ Z\_{5} \\ \hline \\ \begin{pmatrix} -4R\_{d}\mathbf{K}\_{6}(1+(\theta\_{w}-1)Z\_{4})^{2}(\theta\_{w}-1)Z\_{5}\mathbf{Z}\_{5} - \mathbf{K}\_{7}((Z\_{1}\mathbf{Z}\_{5}-2Z\_{2}\mathbf{Z}\_{4}) - \varepsilon(2\mathbf{Z}\_{4} + \frac{\eta}{2}\mathbf{Z}\_{5})) - \\ \hline \\ \begin{pmatrix} 1+\frac{4}{3}R\_{d}\mathbf{K}\_{6}(1+(\theta\_{w}-1)Z\_{4})^{3} \end{pmatrix} \end{pmatrix} \end{pmatrix} \tag{29}$$

with initial conditions (ICs) as follows:

$$\begin{pmatrix} Z\_1(0) \\ Z\_2(0) \\ Z\_2(\infty) \\ Z\_5(0) \\ Z\_4(\infty) \end{pmatrix} = \begin{pmatrix} 0 \\ 1 \\ 0 \\ -\mathcal{K}\_6 \mathfrak{f} (1 - Z\_4(0)) \\ 0 \end{pmatrix} . \tag{30}$$

Similarly,

$$\frac{d}{d\eta} \begin{pmatrix} Z\_1 \\ Z\_2 \\ Z\_3 \\ Z\_4 \\ Z\_5 \end{pmatrix} = \begin{pmatrix} Z\_2 \\ Z\_3 \\ \frac{-\left\{\mathbf{K}\_2(\mathbf{Z}\_1 \mathbf{Z}\_3 - \mathbf{Z}\_2 \mathbf{Z}\_2 - \varepsilon \left(\frac{\mathbf{Z}}{2} \mathbf{Z}\_3 + \mathbf{Z}\_2\right)\right\} - \mathrm{M} \mathbf{K}\_3 \mathbf{Z}\_2 + \mathbf{K}\_4 \mathbf{W} \mathbf{Z}\_4\right\}}{\mathbf{K}\_5} \\\ Z\_5 \\\ \begin{pmatrix} -\mathbf{4} \mathbf{R}\_4 \mathbf{K}\_8 \left(1 + \left(\theta\_w - 1\right) \mathbf{Z}\_4\right)^2 \left(\theta\_w - 1\right) \mathbf{Z}\_5 \mathbf{Z}\_5 - \mathbf{K}\_9 \left(\left(\mathbf{Z}\_1 \mathbf{Z}\_5 - 2 \mathbf{Z}\_2 \mathbf{Z}\_4\right) - \varepsilon \left(2 \mathbf{Z}\_4 + \frac{\eta}{2} \mathbf{Z}\_5\right)\right) - \\\ \frac{\mathbf{K}\_6 \left(\mathbf{A}\_6 \mathbf{Z}\_6 + \theta\_6 \mathbf{Z}\_6\right) - \mathbf{P}\_V \mathbf{K}\_8 \mathbf{Z}\_5 \mathbf{Z}\_3}{\left(1 + \frac{\eta}{4} \mathbf{R}\_4 \mathbf{A}\_8 \left(1 + \left(\theta\_w - 1\right) \mathbf{Z}\_4\right)^3\right)} \end{pmatrix}, \end{pmatrix}, \tag{31}$$

with the corresponding (ICs) as follows:

$$\begin{pmatrix} Z\_1(0) \\ Z\_2(0) \\ Z\_2(\infty) \\ Z\_5(0) \\ Z\_4(\infty) \end{pmatrix} = \begin{pmatrix} 0 \\ 1 \\ 0 \\ -\mathcal{K}\_8 \mathfrak{f}(1 - Z\_4(0)) \\ 0 \end{pmatrix} . \tag{32}$$

The use of an efficient estimation for *F* <sup>00</sup>(0) and *θ* 0 (0) until the boundary restriction is reached addresses these equations. The step size is fixed to ∆*η* = 0.01, which is sufficient to achieve the graphical and the numerical result in tabular form. The range is taken to be *η*max = 10, where the finite value of the dimensional variable *η* for the boundary restrictions is *η*max. The convergence criteria and the accuracy of the outcomes in all cases are up to level 10−<sup>10</sup> .
