**4. Results and Discussion**

The bvp4c function in Matlab is adopted to produce the results of the nonlinear system of ordinary di fferential Equations (8)–(10) together with the boundary conditions (11). The relative error tolerance is set as 10−<sup>10</sup> to gain the results of the numerical outcomes and stability analysis. Table 3 presents that the average central processing unit (CPU) time required for computing each result in Table 4 is approximately 1.4 s, and the dual solutions were obtained by indicating di fferent initial guesses which must satisfy the far-field boundary conditions (11) asymptotically. The numerical results are validated with the numerical results produced by Mahapatra and Sidui [10] and Nawaz and Hayat [57] by collating the values of the shear stress *f* (0) of an axisymmetric (γ = 0) stagnation point flow of a viscous fluid with the exclusion of magnetic field and the unsteady parameter, which is depicted in Table 4. It demonstrates excellent agreemen<sup>t</sup> with the previous literature; hence, the practicality and e ffectiveness of the bvp4c method are verified. The estimated relative error, ε*r* is also measured, and it shows that the calculated values of ε*r* are relatively small between the present and previous results. Figure 2 exhibits the variation of the wall shear stress parameter *f* (0) and *g*(0) towards the di fference value of γ when φ1 = φ2 = *M* = *A* = 0, λ = 0.1, ε = 1.0 and Pr = 1.0 where the dotted lines correspond to the asymptotic behaviour of *f* (0) and *g*(0) which is in excellent agreemen<sup>t</sup> with Mahapatra and Sidui [10] who pursued a standard fourth-order Runge-Kutta integration technique in their study. This justifies the role of the bvp4c numerical technique as a dependable practice, and the present results are valid and correct.


**Table 3.** Computational time to generate the values of *f* (0) as λ varies when φ1 = φ2 = 0, γ= 0, ε= 1, *M* = 0, A = 0, and Pr = 6.2.

The dimensionless velocity profiles *f* (η) and *g* (η) for di fferent values of λ are illustrated in Figures 3 and 4. The figures prove that the profiles of the velocity comply with the far-field boundary conditions of Equation (11) asymptotically. The maximum value of the velocity gradient with the lowest thickness of the momentum boundary layer is preserved for the largest value of λ. Notably, the distance of two adjacent profiles increases remarkably with the increased amount of λ in both the first and second solutions. Figures 5 and 6 portray the variations of *f* (0) and *g*(0) towards ε for di fferent values of λ in hybrid nanofluids with the existence of the magnetic field and the influences of the unsteadiness parameter. An enrichment in the amounts of λ embarking on the augmentation of both λ + γ and λ − γ provide a significant e ffect to the surface shear stresses. The variations of *f* (0) are expected to be higher by the increasing value of λ + γ ye<sup>t</sup> decrease when λ + γ is decreased, and the same trend is expected in *g*(0) with λ − γ. On the other hand, the e ffect of the nanoparticles volume fraction is observed in Figures 7–9, respectively. Surprisingly, the critical values of the various usage of the nanoparticles' volume fraction give no significant e ffect to the trend of the nanofluid (φ1 = 0.02, φ2 = 0.00) and hybrid nanofluid flow (φ1 = 0.02, φ2 = 0.002, 0.04). It is observed that the

reduced skin friction coefficient in *x* and *y* directions, as presented in Figures 7 and 8, and the reduced local Nusselt number in Figure 9 upsurges with the rising in the nanoparticles' volume fraction. This is due to the fact that more kinetic energy is produced with a higher concentration of nanoparticles and thus, enhances the heat transfer of the fluid particles. This finding also corresponds with several present particle-laden direct numerical studies (DNS), which reveal that the roughness components appear to redistribute the energy and, therefore, decrease the overall large-scale near-wall anisotropy of the flow pattern (see Yuan and Piomelli [58] and Ghodke and Apte [59]).

**Table 4.** Results of *f* (0) for specific values of λ when φ1 = 0.0, φ2 = 0.0, γ= 0, ε= 1, *M* = 0, A = 0, and Pr = 6.2.


<sup>∗</sup>ε*r* is the estimate percentage relative error between the present result, *r* and previous result, *s*.

**Figure 2.** Variations of *f* (0) and *g*(0) for different values of γ.

**Figure 3.** Velocity profiles of *f*(η) for different values of λ.

**Figure 4.** Velocity profiles of *g*(η) for different values of λ.

**Figure 5.** Variations of *f* (0) for different values of λ.

**Figure 6.** Variations of *g*(0) for different values of λ.

**Figure 7.** Variations of *f* (0) for different values of φ2.

**Figure 8.** Variations of *g*(0) for different values of φ2.

Figures 10–12 exhibit the influence of the magnetic parameter on *f* (0), *g*(0) and −θ(0) which shows a prominent effect on the fluid flow of the shrinking sheet. The reduced skin friction coefficient in both the *x*− and *y*− directions rise and eventually increase the value of −θ(0) with the escalation of *M* due to the occurrence of the Lorentz force, which acts to retard the fluid flow. The Lorentz force creates resistance to the motion of the fluid particles and then consequently reduces the fluid velocity. The synchronism of the magnetic and electric field that occurred from the formation of the Lorentz force tends to slow down the fluid movement. The boundary layer becomes thinner in both directions as proven in Figures 13 and 14 as *M* increases due to the delayed flow, hence contributing to the increment of *f* (0) and *g*(0). On the other hand, Figure 15 advertises the variations of the non-dimensional temperature profiles for different values of *M* which tend to decrease in the first solution but show a reverse trend in the second solution. The figures clearly reveal that the thicknesses of the thermal boundary layers decrease with the increase of *M* Practically, by restricting the magnetic field intensity, the progression of the thermal boundary layers' thicknesses can be managed and, thus, be able to reduce the temperature profile distributions.

**Figure 9.** Variations of −θ(0) for different values of φ2.

**Figure 10.** Variations of *f* (0) for different values of *M*.

**Figure 11.** Variations of *g*(0) for different values of *M*.

**Figure 12.** Variations of −θ(0) for different values of *M*.

**Figure 13.** Velocity profiles of *f*(η) for different values of *M*.

**Figure 14.** Velocity profiles of *g*(η) for different values of *M*.

The effect of the unsteadiness parameter *A* on *f* (0), *g*(0) and −θ(0) is depicted in Figures 16–18, respectively. Increasing the values of *A* results in a reduction of the skin friction coefficients in both directions, as illustrated in Figures 16 and 17, which consequently decreases the reduced local Nusselt number. The fact that the unsteadiness parameter affects the velocity and temperature profile is proven in Figures 19–21. The velocity of the fluid is in decline along the surface of the sheet due to the increasing value of the shear stress and subsequently shrinks the thickness of the momentum boundary layer nearby the wall, as displayed in Figures 19 and 20, accordingly. The increasing value

of *A* decreases in the temperature profile of the hybrid nanofluid, as shown in Figure 21, which is understandable because the spaces between the molecules become higher in unsteady flow and proportionately decrease the temperature profile and improve the cooling rates of the fluid. Therefore, the unsteadiness parameter should be highlighted and well considered for practical purposes.

**Figure 15.** Temperature profiles of <sup>θ</sup>(η) for different values of *M*.

**Figure 16.** Variations of *f* (0) for different values of *A*.

**Figure 17.** Variations of *g*(0) for different values of *A*.

**Figure 18.** Variations of −θ(0) for different values of *A*.

**Figure 19.** Velocity profiles of *f*(η) for different values of *A*.

**Figure 20.** Velocity profiles of *g*(η) for different values of *A*.

**Figure 21.** Temperature profiles of <sup>θ</sup>(η) for different values of *A*.

Since the dual solutions noticeably exist, as illustrated in Figures 3–21, a stability analysis is performed to discover a significant solution between the first and second solutions. This process is achievable by clarifying the eigenvalue problems in Equations (25)–(27) using bvp4c in Matlab, which produce an infinite set of ω1 < ω2 < ω3 ....... The stability of the flow is reliant on the smallest positive eigenvalue ω1 since there exists an initial decay of disturbances. In contrast, there is an initial growth of perturbations if the smallest eigenvalue ω1 is negative, which signifies that the solution is unstable. As depicted in Table 5, the lowest eigenvalues for the first solution are positive, while the second solution is negative. In conclusion, the first solution is stable and physically reliable, while the second solution is unstable and unreal.


**Table 5.** Smallest eigenvalues ω1 when φ1 = φ2 = 0.02, λ = 1.5, γ = 0.5,*M* = 0.5, A = 0.5 and P*r* = 6.2.
