*2.3. Solution Methodology*

Choosing different variables below

$$y\_1 = \eta, \ y\_2 = f, \ y\_3 = f', \ y\_4 = f'', \ y\_5 = \theta, \ y\_6 = \theta' \tag{14}$$

Equations (7) and (8) can be transformed into a system of first-order differential equations. Then, the Caputo fractional-order derivative is applied to the resultant system to produce a fractional-order system of the following form:

$$\begin{cases} D^{a}\_{\eta}y\_{1} = 1, D^{a}\_{\eta}y\_{2} = y\_{3\prime}, D^{a}\_{\eta}y\_{3} = y\_{4\prime}, D^{a}\_{\eta}y\_{4} = \left(y\_{3}\right)^{2} - y\_{2}y\_{4} + A\left(y\_{3} + \frac{y\_{1}y\_{4}}{2}\right) + My\_{3} \\\ D^{a}\_{\eta}y\_{5} = y\_{6\prime}, D^{a}\_{\eta}y\_{6} = Pr\left[my\_{3}y\_{5} - y\_{2}y\_{6} + A\left(y\_{5} + \frac{y\_{1}y\_{6}}{2}\right) - Ec(y\_{4})^{2}\right] \end{cases} \tag{15}$$

with boundary conditions below:

$$y\_1 = 0, \ y\_2 = \text{S}, \ y\_3 = -1, \ y\_4 = u\_1, \ y\_5 = 1, \ y\_6 = u\_2 \tag{16}$$

Now, the Adams types predictor-corrector method has been applied to get the solution of fractional differential equations. The error of this method is of the order *h*5, where *h* is the grid size.

#### **3. Results and Discussion**

The governing Equation (15) with initial conditions (16) is solved using the Adams-type predictor-corrector method, and dual solutions are found depending on the suction parameter. It is worth mentioning that when an unsteady parameter is equal to 0, our equation of momentum is reduced to an equation obtained in [23], which is the major reference of our work. Furthermore, the results of the coefficient of skin friction in our problem are in good agreement with their work. According to our results, the coefficient of skin friction is equal to *f* --(0) = *<sup>S</sup>*<sup>±</sup> √ *S*2−4+4*M* <sup>2</sup> . For details on the comparison, please refer to Table 1. From the table, we can conclude that only the first solution of our problem is a physical realizable solution, since *f* --(0) > 0, whereas the second solution is unstable because most of the values of *f* --(0) are less than 0. It is worth mentioning that our results of the first solution are approximately equal to the result of a published article [23], which gives us confidence on our calculation (see Table 1). It should be noted that if α < 1, multiple solutions do not exist and do not fulfill the boundary conditions asymptotically.


**Table 1.** Comparison *f* --(0) of present results with [23].

Figure 1 shows the effect of the suction parameter on the velocity profile. It was noticed that the thickness of the velocity boundary layer decreased as suction increased in the first solution. This occurred due to the fact that high suction produced the resistance in the fluid flow, and, as a result, the velocity and thickness of the momentum boundary layer decreased. On the other hand, the suction was proportional to the velocity profile in the second solution.

**Figure 1.** *f*- (η) for increasing values of *S*.

The effect of the magnetic parameter *M* on the velocity is demonstrated in Figure 2. The velocity of the fluid flow decreased as magnetic parameter *M* increased in the first solution, as expected. This was due to the Lorentz or electromagnetic force, which can be defined as "the force which is exerted by a magnetic field on a moving fluid" [39]. We can say this force opposes the transport phenomenon. However, the opposite trend can be seen in the second solution.

**Figure 2.** *f*- (η) for increasing values of *M* and α = 1.

Based on the results shown in Figure 3, there no change was noticed in the first solution when the magnitude of the unsteadiness parameter increased. On the other hand, the velocity layer became thicker initially and then thinner in the second solution, since deaccelerating of the unsteadiness parameter produced more drag force, which caused the thickness of the momentum boundary layer to decrease.

**Figure 3.** *f*- (η) for increasing values of *A* and α = 1.

The effect of the unsteadiness parameter *A* on a dimensionless profile of temperature is depicted in Figure 4. Both thermal boundary layer thicknesses and temperatures decreased initially and then started to increase when the unsteadiness parameter *A* was increased in the second solution. This behavior was expected because the momentum boundary layer declined and, therefore, the temperature increased. However, no difference could be seen in the first solution with the increasing magnitude of the unsteadiness parameter *A*.

**Figure 4.** θ(η) for increasing values of *A* and α = 1.

Figure 5 was drawn for the Prandtl number effect *Pr* on the profile of temperature. We can see that the temperature declined with respect to *Pr* = 0.04 to 6.7 in the first solution. This was because "fluid has relatively lower thermal conductivity for a large value of Prandtl number, which decreases the conduction and thickness of the thermal boundary layer" [40], and, consequently, the temperature reduced. This was because the "Prandtl number *Pr* which is the ratio of momentum diffusion to thermal" [18]. On the other hand, the thermal boundary layer thickness and temperature increase in the range of 0.04 ≤ *Pr* ≤ 3 decreased in the range of 3 ≤ *Pr* ≤ 6.7 in the second solution.

**Figure 5.** θ(η) for increasing values of *Pr* and α = 1.

Figure 6 indicates the temperature increase when the Ecker number was increased in the first solution. This was due to fact that an expansion in dissipation enhanced the flow of thermal conductivity, which extended the temperature and thermal boundary layer thickness. On the other hand, the Eckert number was inversely proportional to the temperature and thickness of the thermal boundary layer in the second solution.

**Figure 6.** θ(η) for increasing values of *Ec* and α = 1.

The graph of the coefficient of skin friction for several values of *S* and different values of *A* is illustrated in Figure 7. It was observed that the skin friction coefficient increased (decreased) when suction was increased (reduced) in the first (second) solution. However, it decreased with the decreasing of the unsteady parameter *A*. Physically, resistance occurred due to increments in the suction parameter in the stable solution, while the opposite trend was seen in the unstable solution.

**Figure 7.** Coefficient of skin friction for different values of *A*.

Figure 8 shows the effect of α on the profile of the temperature. It was noticed that multiple solutions were difficult to be obtained when α < 1. As α increased, the temperature of the fluid increased in the first solution and decreased in the second solution.

**Figure 8.** θ(η) for increasing values of α.

Figure 9 demonstrates the effect of α on the velocity profile. In both solutions, the velocity of the fluid decreased when α increased.

**Figure 9.** *f*- (η) for increasing values of α.
