*4.2. Three-Stage Lobatto III-A Formula*

Three-stage Lobatto III-A formula is built in BVP4C function with aid of *C* <sup>1</sup> piece-wise cubic polynomial in the finite difference code. According to Lund et al., [35] and Raza et al., [36], "this collocation polynomial and formula offers a *C* 1 continuous solution in which mesh error control and selection are created on the residual of the continuous solution. The tolerance of relative error is fixed 10−<sup>5</sup> for the current problem. The suitable mesh determination is created and returned in the field sol.x. The bvp4c returns solution, called as sol.y., as a construction. In any case, values of the solution are gotten from the array named sol.y relating to the field sol.x". In addition, Figure 2 explains the algorithm of the method for stability analysis of the solutions.

*Crystals* **2020**, *10*, x FOR PEER REVIEW 7 of 14

**Figure 2.** Description of three-stage Lobatto III-A formula. **Figure 2.**Description of three-stage Lobatto III-A formula.

#### **5. Result and Discussion 5. Result and Discussion**

Aurangzaib et al., [29] solved these Equations (11)–(14) and found dual solutions, which was the main contribution of authors. Aurangzaib et al., [29] have given some strong statements during the investigation of critical points. Before going to detail, we check the accuracy of our method (shooting method) by comparing our results with previously published literature in Table 1 and found excellent agreement with them. Moreover, we also qualitatively compared our results with Aurangzaib et al., [29] in Figure 3 and found in the good agreement. This gives trust in our numerical calculation and urges us to additionally contemplate this problem. Aurangzaib et al., [29] found dual solutions and stated in the result and discussion section that "the present study shows that for *K* = 0.1, i.e., for micropolar fluid, the similarity solutions exist when S ≥ 2.3231 and no similarity solution exists for S < 2.3231". Firstly, we would like to clarify that there exist triple solutions not dual solutions; secondly, there is a range of multiple solutions and single solution. From Figure 4, a conclusion can be made that there exists multiple similarity solutions when S ≥ 2.3224 and only a single similarity solution exists for S < 2.3224 when *K* = 0.1. However, the range of multiple similarity solution is 2.3769 ≤ S and there also exists only a single solution when S < 2.3769 when *K* = 0.2 (see Figure 4). Furthermore, skin friction increases as suction is increased in the third solution. It is worth mentioning here that there is no range of no solution. This is one of the big reasons that insist us to reconsider and re-examine the whole problem by reproducing all of the results because there are triple solutions in order to provide true knowledge to the readers and researchers. Aurangzaib et al., [29] solved these Equations (11)–(14) and found dual solutions, which was the main contribution of authors. Aurangzaib et al., [29] have given some strong statements during the investigation of critical points. Before going to detail, we check the accuracy of our method (shooting method) by comparing our results with previously published literature in Table 1 and found excellent agreement with them. Moreover, we also qualitatively compared our results with Aurangzaib et al., [29] in Figure 3 and found in the good agreement. This gives trust in our numerical calculation and urges us to additionally contemplate this problem. Aurangzaib et al., [29] found dual solutions and stated in the result and discussion section that "the present study shows that for *K* = 0.1, i.e., for micropolar fluid, the similarity solutions exist when S ≥ 2.3231 and no similarity solution exists for S < 2.3231". Firstly, we would like to clarify that there exist triple solutions not dual solutions; secondly, there is a range of multiple solutions and single solution. From Figure 4, a conclusion can be made that there exists multiple similarity solutions when S ≥ 2.3224 and only a single similarity solution exists for S < 2.3224 when *K* = 0.1. However, the range of multiple similarity solution is 2.3769 ≤ S and there also exists only a single solution when S < 2.3769 when *K* = 0.2 (see Figure 4). Furthermore, skin friction increases as suction is increased in the third solution. It is worth mentioning here that there is no range of no solution. This is one of the big reasons that insist us to reconsider and re-examine the whole problem by reproducing all of the results because there are triple solutions in order to provide true knowledge to the readers and researchers.

**Table 1.** The comparison of values of heat transfer rate for different values of *Pr*.


**Figure 3.** Comparison between upper graph: existing results obtained by Aurangzaib et al., [29] and lower graph: present results. **Figure 3.** Comparison between upper graph: existing results obtained by Aurangzaib et al., [29] and lower graph: present results. **Figure 3.** Comparison between upper graph: existing results obtained by Aurangzaib et al., [29] and lower graph: present results.

**Figure 4.** Variation of ݂′′(0) for different fixed values of *K* with suction parameter *S*. **Figure 4.** Variation of ݂′′(0) for different fixed values of *K* with suction parameter *S*. **Figure 4.** Variation of *f* 00(0) for different fixed values of *K* with suction parameter *S*.

Figure 5 illustrates the effect of micropolar parameter on the *h* 0 (0). The effect of local couple stress enhanced as the suction increases in the first and third solutions because increasing suction creates additional resistance in the flowing fluid inside the boundary layer. However, increments in the material parameter produce more coupling of stress. Figure 6 shows the nature of the heat transfer rate for various values of the suction. It has been examined that the heat transfer rate increases in the first and second solutions for the higher values of the suction parameter, while the third solution shows opposite compliance. Figure 5 illustrates the effect of micropolar parameter on the ℎ′(0). The effect of local couple stress enhanced as the suction increases in the first and third solutions because increasing suction creates additional resistance in the flowing fluid inside the boundary layer. However, increments in the material parameter produce more coupling of stress. Figure 6 shows the nature of the heat transfer rate for various values of the suction. It has been examined that the heat transfer rate increases in the first and second solutions for the higher values of the suction parameter, while the third solution shows opposite compliance. Figure 5 illustrates the effect of micropolar parameter on the ℎ′(0). The effect of local couple stress enhanced as the suction increases in the first and third solutions because increasing suction creates additional resistance in the flowing fluid inside the boundary layer. However, increments in the material parameter produce more coupling of stress. Figure 6 shows the nature of the heat transfer rate for various values of the suction. It has been examined that the heat transfer rate increases in the first and second solutions for the higher values of the suction parameter, while the third solution shows opposite compliance.

*Crystals* **2020**, *10*, x FOR PEER REVIEW 9 of 14

*Crystals* **2020**, *10*, x FOR PEER REVIEW 9 of 14

**Table 1.** The comparison of values of heat transfer rate for different values of *Pr*.

*Pr M* **Ishak [33] Pramanik [37] Raju et al. [38] Present Results** 0 0.9548 0.9547 0.954734 0.954955 0 1.4715 1.4714 1.471426 1.471421 0 1.8691 1.8691 1.869134 1.869044

**Table 1.** The comparison of values of heat transfer rate for different values of *Pr*.

*Pr M* **Ishak [33] Pramanik [37] Raju et al. [38] Present Results** 1 0 0.9548 0.9547 0.954734 0.954955 2 0 1.4715 1.4714 1.471426 1.471421

10 0 3.6603 3.6603 3.660312 3.660354

**Figure 5.** Variation of ℎ′(0) for different fixed values of *K* with suction parameter *S*. **Figure 5.** Variation of *h* 0 (0) for different fixed values of *K* with suction parameter *S*. **Figure 5.** Variation of ℎ′(0) for different fixed values of *K* with suction parameter *S*.

**Figure 6.** Variation of −ߠ′)0 (for different fixed values of *K* with suction parameter *S*. **Figure 6.** Variation of −θ 0 (0) for different fixed values of *K* with suction parameter *S*.

**Figure 6.** Variation of −ߠ′)0 (for different fixed values of *K* with suction parameter *S*. Finally, we plot Figures 7–10 to show the existence of triple solutions of velocity, microrotation, and temperature profiles for different values of material parameter *K*. In Figure 7, the dual nature of behavior has been noticed in the first solution. Velocity profile increases in the second solution when *K* is increased; the physical material parameter reduces the effect of drag force due to that thickness of the momentum boundary layer enhanced. On the other hand, the opposite trend has been observed in the third solution. Dual behavior can be noticed microrotation profile in Figure 8

for all solutions. The thickness of thermal boundary layer increases in the first and second solutions as material parameter *K* is increased, as in Figure 9, since the non-Newtonian parameter produces more viscosity, decreases the velocity of profiles, and forces fluid flow to stay on the hotter surface, as a result temperature of fluid increases and the boundary layer becomes thicker. However, the opposite behavior is noticed for the third solution. The Prandtl effect on the temperature distribution is depicted in Figure 10. It is observed that the temperature of fluid diminishes for the higher values of the Prandtl number for all solutions. Physically, it can be explained, as the Prandtl number (*Pr* = µ*cp k* ) has an inverse relationship with the thermal conductivity and, consequently, diminishes the thickness of the thermal boundary layer. The thickness of thermal boundary layer increases in the first and second solutions as material parameter *K* is increased, as in Figure 9, since the non-Newtonian parameter produces more viscosity, decreases the velocity of profiles, and forces fluid flow to stay on the hotter surface, as a result temperature of fluid increases and the boundary layer becomes thicker. However, the opposite behavior is noticed for the third solution. The Prandtl effect on the temperature distribution is depicted in Figure 10. It is observed that the temperature of fluid diminishes for the higher values of the Prandtl number for all solutions. Physically, it can be explained, as the Prandtl number (ܲݎ= ఓ ) has an inverse relationship with the thermal conductivity and, consequently, diminishes the thickness of the thermal boundary layer.

in the third solution. Dual behavior can be noticed microrotation profile in Figure 8 for all solutions.

*Crystals* **2020**, *10*, x FOR PEER REVIEW 10 of 14

Finally, we plot Figures 7–10 to show the existence of triple solutions of velocity, microrotation, and temperature profiles for different values of material parameter *K*. In Figure 7, the dual nature of behavior has been noticed in the first solution. Velocity profile increases in the second solution when

**Figure 7.** Effect of *K* on ݂′(ߟ (when *S = 2.5*. **Figure 7.** Effect of *K* on *f* 0 (η) when *S* = *2.5*. *Crystals* **2020**, *10*, x FOR PEER REVIEW 11 of 14

**Figure 8.** Effect of *K* on ℎ(ߟ (when *S = 2.5*. **Figure 8.** Effect of *K* on *h*(η) when *S* = *2.5*.

**Figure 9.** Effect of *K* on ߠ)ߟ (when *S = 2.5.*

**Figure 8.** Effect of *K* on ℎ(ߟ (when *S = 2.5*.

**Figure 10.** Effect of *Pr* on ߠ)ߟ(. **Figure 10.** Effect of *Pr* on θ(η).

Table 2 shows the smallest eigenvalues ߝ for the selected values of *S* and *K*. The positive smallest eigenvalue makes the initial disturbance decay and, in this way, the flow becomes stable. Conversely, the negative smallest eigenvalue outcomes in an initial growth of disturbance, in this manner, the flow is unstable. It is seen from Table 2 that ߝ is negative for the second and the third solutions, while positive for the first solution. Thus, the second and the third solution are not stable, and the first solution is stable. From this discussion, it can be concluded that the first solution of Aurangzaib et al., [29] is not stable and not physically realizable; therefore, in this stage, it could be said that the first solution is actually the second or the third solution. Table 2 shows the smallest eigenvalues ε for the selected values of *S* and *K*. The positive smallest eigenvalue makes the initial disturbance decay and, in this way, the flow becomes stable. Conversely, the negative smallest eigenvalue outcomes in an initial growth of disturbance, in this manner, the flow is unstable. It is seen from Table 2 that ε is negative for the second and the third solutions, while positive for the first solution. Thus, the second and the third solution are not stable, and the first solution is stable. From this discussion, it can be concluded that the first solution of Aurangzaib et al., [29] is not stable and not physically realizable; therefore, in this stage, it could be said that the first solution is actually the second or the third solution.

The micropolar fluid over the shrinking surface has been considered. The system of governing equations has been transformed into the system of ODEs by using appropriate exponential similarity transformation. The system of ODEs is reduced to IVPs by employing the shooting method before solving IVPs by the Runge Kutta method. The pointwise conclusions of this study are given below:

2. There are ranges of multiple solutions and no solutions that depend upon the suction

3. According to stability analysis, the first solution is stable, which can be experimentally seen.

5. The thickness of thermal boundary layer increases in the first and the second solutions as

st solution 2

0.1 2.3224 1.28061 0 0 - 2.4 1.0662 −0.06382 −0.13406 0.2 2.3769 1.36201 0 0 - 2.4 1.1364 −0.10482 −0.17482

ࢿ

nd solution 3

rd solution

*K S*

4. The results of Aurangzaib et al., [29] are unstable.

6. Increments in the material parameter produce more couple stress.

material parameter *K* is increased.

**6. Conclusion**

1. Triple solutions appear.

parameter.

1


**Table 2.** Smallest eigen values ε at several values of *S* and *K* when *Pr* = 1.
