**4. Numerical Procedure**

The nonlinear coupled ODEs (9) and (10) with boundary constraint (11) through the bvp4c by converting the leading ODEs to an initial value problem (IVP). In this method, it is further helpful to provide a fixed value to η → <sup>∞</sup>, say η<sup>∞</sup>. The above-mentioned higher order equations are converted into a first-order system as follows:

$$f' = p\_\prime \tag{15}$$

$$p' = q\_{\prime} \tag{16}$$

$$q' = \frac{\left(1 + \left(\mathcal{W}cf'\right)^{1-n}\right)^2}{\varepsilon\_1 \left(1 + n\left(\mathcal{W}cf'\right)^{1-n}\right)^{1-n}} \Big| p^2 - fq - 1 + K(1-p) - B\left(1 - p^2\right) - \lambda\theta \Big|\_{\prime} \tag{17}$$

$$
\theta' = \underline{z},
\tag{18}
$$

$$z'=p\theta-fz,\tag{19}$$

with

$$f(0) = 0, \; p(0) = \gamma\_1 q(0), \; \theta(0) = 1 + \gamma\_2 z(0). \tag{20}$$

Numerically grip the system of Equations (15)–(20) as an IVP, we require that the values for *q*(0) and θ(0) are needed, however these values are not mentioned. The initial estimated values for *q*(0) and θ(0) are conjectured and bvp4c is pertained on MATLAB software to achieve accurate results. It is also noted that the multiple solutions are attained by setting different guesses. After that, the considered values of <sup>θ</sup>(η) and *f*(η) at (η∞ = 8) are evaluated with the boundary conditions <sup>θ</sup>(η∞) = 0 and *f*(η∞) = 1, in which the predictable values of *q*(0) and θ(0) are prescribed by the Secant method to achieve a better guess for the solutions. The step size is considered as Δη = 0.01. The procedure is iteratively repeated until required solutions with an acceptable level of accuracy (i.e., up to <sup>10</sup>−5) to fulfill the criterion of convergence.
