**4. Numerical Solutions**

To solve the boundary value problems (11)–(13) with boundary conditions given by (15), we have adopted a built-in function called bvp4c from Matlab package. Further, to access the precision of current algorithm, the current results of skin friction coefficient *f* (0) are compared with previously reported solutions of Ishak et al. [10], who used Keller-box method in their work, Mahapatra and Nandy [65], who pursued the shooting method for their computation and Zainal et al. [59] which employed the bvp4c solver. These comparative solutions are revealed in Table 3 for selected values of shrinking parameter (*λ* < <sup>0</sup>). It can be point out from these tables that there is good agreemen<sup>t</sup> with these methods (error is relatively small), thereby confirming the consistency of the approach used. Furthermore, this validates the present model and proves the accuracy of the bvp4c solver in solving a boundary layer problem as the present results able to withstand the Keller-box method and shooting method which have been employed by Ishak et al. [10] and Mahapatra and Nandy [65]. In this part, the results of local skin friction *Cf*Re1/2 *x* , Nusselt number *Nux*Re−1/2 *x* , velocity profile *f* (*η*), microrotation profile *h*(*η*) as well as temperature profile *<sup>θ</sup>*(*η*) are illustrated graphically to explore the influence of some governing parameters such as Cu nanoparticle volume fraction (*ϕ*2), unsteady parameter (*A*), material parameter (*K*) and radiation parameter (*Rd*).



'[ ]' Second solution.

The effect of Cu nanoparticle volume fraction *ϕ*2 against stretching or shrinking parameter *λ* on the local skin friction *Cf*Re1/2 *x* and Nusselt number *Nux*Re−1/2 *x* as given in Equation (17) are shown in Figure 2a,b. It is apparent from these figures that for shrinking parameter (*λ* < <sup>0</sup>), the occurrence of dual solutions is noted. However, it is remarked that no solution exists when *λ* < *λ<sup>c</sup>*, which indicates that the boundary layer is detach from the surface and the principle of boundary layer theory are no longer valid. Moreover, *λc* is the critical point that connected the first and second solutions. In addition, a unique solution is noticed when the sheet is stretched (*λ* > <sup>0</sup>). It is clear that for *ϕ*2 = 0, the problem reduces to the micropolar nanofluid. From these figures, it is discovered that upsurge in Cu nanoparticle volume fraction *ϕ*2 enhances the local skin friction *Cf*Re1/2 *x* and local Nusselt number *Nux*Re−1/2 *x* for all domains of stretching and shrinking parameter *λ* in first solution but a small change is noted for the second solution. This finding proves that the increment of Cu nanoparticle volume fraction *ϕ*2 can improve the heat transfer efficiency. This also implies that hybrid nanofluid provides a better heat performance than nanofluid. Furthermore, the enhancement of Cu nanoparticle volume fraction *ϕ*2 on local skin friction *Cf*Re1/2 *x* and Nusselt number *Nux*Re−1/2 *x* fasten the detachment of boundary layer flow. Figure 3a–c exemplify the impact of Cu nanoparticle volume fraction parameter *ϕ*2 on the velocity profile *f* (*η*), microrotation profile *h*(*η*) and temperature profile *<sup>θ</sup>*(*η*) for shrinking sheet (*λ* = <sup>−</sup>1.25). It reveals that augmentation of Cu nanoparticle volume fraction *ϕ*2 depreciates the momentum and microrotation boundary layer thickness in first and second solutions. Meanwhile, the thermal boundary layer thickness increases as Cu nanoparticle volume fraction *ϕ*2 increases for the first solution, however a contrary observation is noted for the second solution. Furthermore, one can see that the boundary layer thickness of the first solution was slimmer than second solution. In addition, all the profiles published are asymptotically satisfied the boundary conditions (15) and eventually supported the findings obtained in Figure 2a,b.

**Figure 2.** (**a**) *Cf* Re1/2 *x* ; (**b**) *Nux*Re−1/2 *x* with *λ* for various *ϕ*2.

**Figure 3.** (**a**) *f* (*η*); (**b**) *h*(*η*); (**c**) *<sup>θ</sup>*(*η*) for various *ϕ*2.

Figure 4a,b show the impact of unsteady parameter *A* = 0, 0.1, 0.2 on the local skin friction *Cf*Re1/2 *x* and Nusselt number *Nux*Re−1/2 *x* towards stretching/shrinking parameter *λ*. The occurrence of unsteady parameter *A* consequently elevates the local skin friction *Cf*Re1/2 *x* and Nusselt number *Nux*Re−1/2 *x* . It must be noted that the flow corresponds to the steady micropolar flow when *A* = 0 and it is numerically observed that dual solution exists as −1.24641 < *λ* < −1 and a unique solution exists when *λ* ≥ −1. However, the physical character of flow changes when the flow becomes unsteady. For instance, as the unsteadiness parameter increase, i.e., *A* = 0.1 and *A* = 0.2, the range of similarity solutions to exist also increases where dual solutions is observed in the range of −1.31162 < *λ* < −1 and −1.38029 < *λ* < −1, respectively. The unique solution however only exists as *λ* ≥ −1 for both case of *A* and concurrently no solution is noticed when *λ* < *λ<sup>c</sup>*. In short, a raise in unsteady parameter *A* act in postponing the boundary layer detachment. On the other side, Figure 5a–c portray the discrepancy of velocity *f* (*η*), microrotation *h*(*η*) and temperature *<sup>θ</sup>*(*η*) profiles when unsteady parameter *A* fluctuates from 0 to 0.15. For the first solution, the diminished of momentum boundary layer thickness is observed when the unsteady parameter *A* increases as shown in Figure 5a, but a reverse observation is remarked for the second solution. Furthermore, the microrotation boundary layer thickness diminishes with an increment of unsteady parameter *A* values in the first solution except near the sheet, while an opposing trend is remarked for the second solution. It is also interesting to observe from these figures that for both solutions, the thermal boundary layer thickness increase with an upsurge of unsteady parameter *A*.

**Figure 4.** (**a**) *Cf* Re1/2 *x* ; (**b**) *Nux*Re−1/2 *x* with *λ* for various *A*.

Figure 6a,b and Figure 7a–c are portray to discuss the effect of material parameter *K* on the local skin friction *Cf*Re1/2 *x* , local Nusselt number *Nux*Re−1/2 *x* , velocity profile *f* (*η*), microrotation profile *h*(*η*) and temperature profile *<sup>θ</sup>*(*η*) for Cu-Al2O3/water. It is obvious that existence of material parameter (*K* = 1, 2) give rises to the local skin friction *Cf*Re1/2 *x* if compared to the absence of material parameter (*K* = <sup>0</sup>), i.e., no vortex viscosity. However, different results are observed for the local Nusselt number *Nux*Re−1/2 *x* where the nonexistence of micropolar fluid (*K* = 0) cause an enhancement in comparison with the existence of material parameter (*K* = 1, <sup>2</sup>). This phenomenon reveals the fact that upsurge value of material parameter gives rise on the vortex viscosity in the fluid flow which consequently enhance the skin friction at the wall and decrease the rate of heat transfer at the wall. Additionally, we observed that an upsurge values of material parameter prompt the domain of similarity solutions to exist become narrow. For instance, the similarity solutions in the nonexistence of material parameter are noted in the range of −1.31178 ≤ *λ* ≤ 1, while in the existence of material parameter (*K* = 1, <sup>2</sup>), the range of solutions are observed to be −1.31164 ≤ *λ* ≤ −1 and −1.31162 ≤ *λ* ≤ −1. Furthermore, an upsurge values of material parameter *K* causing the thickness of the momentum and thermal boundary layer to increase in first and second solutions. We can see from Figure 7b that the microrotation boundary layer thickness for first and second solutions near the

sheet decreases when this parameter rises while the contrary trend is observed for the large *η*.

**Figure 6.** (**a**) *Cf* Re1/2 *x* ; (**b**) *Nux*Re−1/2 *x* with *λ* for various *K*.

**Figure 7.** (**a**) *f* (*η*); (**b**) *h*(*η*); (**c**) *<sup>θ</sup>*(*η*) for various *K*.

The influence of radiation parameter *Rd* on the local Nusselt number *Nux*Re−1/2 *x* and temperature profile *<sup>θ</sup>*(*η*) are portrayed in Figure 8a,b, respectively. We have noticed that the radiation parameter *Rd* has no control on the flow field, which is evident from Equations (11)–(13). The local Nusselt number is discovered to increase with an upsurge values in radiation parameter *Rd*. Furthermore, the thermal boundary layer thickness become thicker when this parameter raises in the first solution and it became slimmer in the second solution. As can be seen from Figure 8b, this parameter however does not give effect on the range of similarity solutions to occur, i.e., the critical value for stretching/shrinking parameter *λc* is the same for all values of radiation parameter *Rd*. It is noteworthy that the distribution of temperature in the fluid is significantly affected by the radiation parameter *Rd*. Physically, this is due to the fact that the heat is produced due to the radiation process and therefore, enhances the fluid temperature.

The boundary value problems (11)–(13) together with boundary conditions (15) observe the occurrence of non–unique solutions for some governing parameters. The phenomenon of non-unique solutions namely first and second solutions are proved and portrayed as in Figures 2–8. Accordingly, an investigation on the stability analysis has been executed in this present work so that we can identified the most stable solutions. Therefore, the linearized Equations (25)–(27) along with conditions (28) have been solved with the aid of the bvp4 function in MATLAB numerically. The smallest eigenvalues *γ* on the selected parameter *A* and *λ* from Figure 4a,b are listed in Table 4. This table shown that the second solution displays the negative values of *γ*, whereas the first solution demonstrates the positive values of *γ*. The smallest eigenvalues *γ* against *λ* have been plotted in Figure 9. This figure eventually supported the findings obtained from Table 4 except when the values of stretching/shrinking parameter *λ* approach its critical value where we observed unstable solutions for both solutions. In reference from the previous literature, we can deduce that the first solution is stable, on the contrary, the second solution is unstable. It worth noting that this analysis is important in identifying the stable solution when there is exist non-unique solutions so that the flow behavior can be predict accurately.


**Table 4.** Smallest eigenvalues *γ* for selected *A* and *λ* when *ϕ*1 = *ϕ*2 = 0.01, *Rd* = 0.1, *K* = 1 and *n* = 0.5.

**Figure 9.** Smallest eigenvalues *γ* against *λ* when *ϕ*1 = *ϕ*2 = 0.01, *A* = *Rd* = 0.1, *K* = 1, *n* = 0.5.
