**4. Numerical Method**

Following [47,48], Equations (23) and (24), subject to (25), will be solved using the local similarity method, where the first derivatives with respect to τ are neglected and the Equations (23) and (24) with boundary conditions (25) can be re-written as

$$2\nabla\_1 f'''' + 2\eta f'' - 4\pi \left( f'^2 - f f'' + \frac{\sigma\_{ff}}{\sigma\_f} \frac{Ha^2}{1 - \phi + \phi(\rho\_s/\rho\_f)} f' \right) = 0 \tag{28}$$

$$\left(\frac{\Xi\_2}{\Pr}\middle|\frac{k\_{ff}}{k\_f} + \frac{4}{3}\mathrm{R}d\right)\theta^{\prime\prime} + 2\eta\theta^{\prime} + 4\tau(f\theta^{\prime} - f^{\prime}\theta)^{\prime} = 0\tag{29}$$

The boundary conditions (25) remain the same. These ordinary di fferential equations with the boundary conditions (25) can be solved numerically by applying the Runge–Kutta–Fehlberg method (RKF7 45). Following [47,48], for the local non-similarity solution, now we hold all the terms by assuming the new auxiliary functions *<sup>F</sup>*(<sup>τ</sup>, η), and <sup>Θ</sup>(<sup>τ</sup>, η), which are defined by

$$F = \frac{\partial f}{\partial \tau}, \Theta = \frac{\partial \theta}{\partial \tau} \tag{30}$$

Thus, Equations (23) and (24) can be expressed as

$$2\pi f'''' + 2\eta f'' - 4\pi \left(f'^2 - ff'' + \frac{\sigma\_{ff}}{\sigma\_f} \frac{Ha^2}{1 - \phi + \phi(\rho\_s/\rho\_f)} f'\right) - 4\pi F' = 0\tag{31}$$

$$2\frac{\Xi\_2}{\text{Pr}} \left(\frac{k\_{ff}}{k\_f} + \frac{4}{3}\text{R}d\right) \theta'' + 2\eta \theta' + 4\tau (f\theta' - f'\theta) - 4\tau \Theta = 0\tag{32}$$

subject to the same condition in (25). The new ODEs (31)–(32), subject to (25) represent a local non-similarity model for the problem under consideration. Equations (31) and (32) and the boundary conditions (25) are now di fferentiated w.r.t. τ, simplified and the derivatives w.r.t. τ are neglected again. These equations represent a local similarity model and can be expressed as;

$$\begin{cases} \Xi\_1 \Gamma'' + 2\eta F' + 4ff' - f'^2 + 4\tau F f' - f'^2 + 4\tau f F' - 2f' F' - \left(\frac{\sigma\_{ff}}{\sigma\_f} \frac{H\alpha^2}{1 - \phi + \phi(\rho\_s/\rho\_f)} F'\right) \\ -4F' = 0 \\ \frac{\Xi\_2}{\Gamma\tau} \left(\frac{k\_{ff}}{k\_f} + \frac{4}{3}\mathcal{R}d\right) \Theta'' + 2\eta \Theta' + 4(f\theta' - f'\theta) + 4\tau F \theta' + f\Theta' - \Theta f' - \theta F' - 4\Theta' = 0 \end{cases} \tag{33}$$

$$\begin{cases} F(\tau,0) = 0, F'(\tau,0) = -\frac{\delta}{2\tau^{3/2}(1-\chi)^{25}} F''(\tau,0) + \frac{\delta/\sqrt{\tau}}{(1-\chi)^{25}} F''(\tau,0), F(\tau,\infty) = 0, \\ \frac{k\_f}{k\_f} \Theta'(\tau,0) = -Bi[\frac{[1-\theta(\tau,0)]}{2\sqrt{\tau}} - \sqrt{\tau}\Theta(\tau,0)], \Theta(\tau,\infty) = 0 \end{cases} \tag{34}$$

The ODEs (31)–(33) subject to (25) and (34) were solved numerically by employing the Runge–Kutta–Fehlberg technique (RKF45) using MAPLE-19 software (MAPLE 2019.0, Maplesoft, Waterloo, ON, Canada). This method is generally known as one of the most excellent methods available for obtaining the solutions of nonlinear di fferential equations and provides more accurate results. The step size was selected. For the similarity variable ηmax, Equations (25) and (34), were replaced as

$$f'(5) = 0, \theta(5) = 0,\\ F'(5) = 0, \Theta(5) = 0 \tag{35}$$

The selection of ηmax = 5 guarantees that all numerical solutions approached the asymptotic values properly.
