**Multiple Solutions of Nonlinear Boundary Value Problems of Fractional Order: A New Analytic Iterative Technique**

**Omar Abu Arqub 1 , Ahmad El-Ajou 1 , Zeyad Al Zhour 2 and Shaher Momani 3,4,\***


*Received: 10 November 2013; in revised form: 11 December 2013 / Accepted: 11 December 2013 / Published: 9 January 2014* 

**Abstract:** The purpose of this paper is to present a new kind of analytical method, the so-called residual power series, to predict and represent the multiplicity of solutions to nonlinear boundary value problems of fractional order. The present method is capable of calculating all branches of solutions simultaneously, even if these multiple solutions are very close and thus rather difficult to distinguish even by numerical techniques. To verify the computational efficiency of the designed proposed technique, two nonlinear models are performed, one of them arises in mixed convection flows and the other one arises in heat transfer, which both admits multiple solutions. Graphical results and tabulate data are presented and discussed quantitatively to illustrate the multiple solutions. The results reveal that the method is very effective, straightforward, and powerful for formulating these multiple solutions.

**Keywords:** multiple solutions; fractional differential equations; residual power series

**AMS Classification:** 32A05; 41A58; 26A33

#### **1. Introduction**

Multiple or dual solutions of nonlinear boundary value problems (BVPs) of fractional order are an interesting subject in the area of mathematics, physics, and engineering. In fact, it is more consequential not to lose any solution of nonlinear BVPs of fractional order due to their wide application in scientific and engineering research. Based on this important fact, the present paper is going to present an analytical method, the so-called residual power series (RPS), that enables us to predict the multiplicity of solutions which nonlinear BVP of fractional order admits and furthermore to calculate the multiple solutions analytically at the same time.

BVPs of fractional order have received considerable attention in the recent years due to their wide applications in the areas of physics and engineering. Many important phenomena in electromagnetics, acoustics, viscoelasticity, electrochemistry, and material science are well described by fractional BVP [1–4]. It is well known that the fractional order differential and integral operators are non-local operators. This is main reason why differential operators of fractional order provide an excellent instrument for description of memory and hereditary properties of various physical and engineering processes. For example, half-order derivatives and integrals proved to be more useful for the formulation of certain electrochemical problems than the classical models [5–9]. Indeed, for example, applying fractional calculus theory to entropy theory has become a significant part and a hotspot research domain [10–19], since the fractional entropy could be used in the formulation of algorithms for image segmentation where traditional Shannon entropy has presented limitations [13] and in the analysis of anomalous diffusion processes and fractional diffusion equations [14–19].

In general, most BVPs of fractional order do not have exact solutions. Particularly, there is no known method for solving these types of equations in closed form solution. As a result, numerical and analytical techniques have been used to study such problems. The reader is referred to [20–27] in order to know more details about the fractional BVPs, including their history and kinds, their existence and uniqueness of solution, their applications and methods of solutions, *etc*.

Series expansions are a very important aid in numerical calculations, especially for quick estimates made in hand calculations, for example, in evaluating functions, integrals, or derivatives. Solutions to differential equations can often be expressed in terms of series expansions. Since, the advent of computers, it has, however, become more common to treat differential equations directly, using different approximation method instead of series expansions, but in connection with the development of automatic methods for formula manipulation, one can anticipate renewed interest in series methods. These methods have some advantages, especially in multidimensional and multiple solutions for BVPs of fractional order.

The RPS method was developed by the first author [28] as an efficient method for determining values of coefficients of the power series solution for first and the second-order fuzzy differential equations. It has been successfully applied in the numerical solution of the generalized Lane-Emden equation, which is a highly nonlinear singular differential equation [29] and in the numerical solution of higher-order regular differential equations [30]. The RPS method is an effective and easy to construct power series solution for strongly linear and nonlinear equations without linearization, perturbation, or discretization [28–30]. Different from the classical power series method, the RPS method does not need to compare the coefficients of the corresponding terms and a recursion relation is not required. This method computes the coefficients of the power series by a chain of equations of one or more variables. In fact, the RPS method is an alternative procedure for obtaining analytic solutions for BVPs of fractional order that admits multiple solutions. By using the residual error concept, we get a series solution, in practice a truncated series solution. Moreover, the multiple solutions and all their fractional derivatives are applicable for each arbitrary point in the given interval. On the other aspect as well, the RPS method does not require any conversion while switching from the low-order to the higher-order; as a result the method can be applied directly to a given problem by choosing an appropriate initial guess approximation.

In the present paper, the RPS method will investigate how to construct new algorithms for predicting and finding multiple solutions for those nonlinear BVPs of fractional order that admit multiple solutions. Furthermore, we will adapt a new generalization of Taylor's series formula that involves Caputo fractional derivatives in order to apply the RPS method.

The results dealing with multiple solutions of BVPs of fractional order are relatively scarce. Recently, many authors have discussed the multiple solutions to some problems using some of the well-known methods. However, the reader is referred to [31–35] in order to know more details about these methods, including their types and history, their motivation for use, their characteristics, and their applications. On the other hand, the numerical solvability of other version of differential equations and other related equations can be found in [36–43] and references therein.

The outline of the paper is as follows: in the next section, we utilize some necessary definitions and results from fractional calculus theory. In Section 3, the general form of generalized Taylor's formula is mentioned and proved. In Section 4, basic idea of the RPS method is presented in order to construct and predict multiple solutions for BVPs of fractional order. In Section 5, two nonlinear models are performed in order to illustrate the capability and simplicity of proposed method. Finally, conclusions are presented in Section 6.

#### **2. Review of Fractional Calculus Theory**

In this section, we present some necessary definitions and essentials results from fractional calculus theory. There are various definitions of fractional integration and differentiation, such as Grunwald-Letnikov's definition and Riemann-Liouville's definition [5,6,8]. The Riemann-Liouville derivative has certain disadvantages when trying to model real-world phenomena with fractional differential equations (FDEs). Therefore, we shall introduce a modified fractional differential operator "j F proposed by Caputo in his work on the theory of viscoelasticity [4]. Throughout this paper,<k the set of integer numbers, l the set of real numbers, and m is the Gamma function.

**Definition 2.1:** A real function \* 5 5 is said to be in the space n\$% ol if there exists a real number p% such that \* 5  5q\* 5, where \* 5 o n) r, and it is said to be in the space n\$ : if \* : 5 o n\$; ok.

**Definition 2.2:** The Riemann-Liouville fractional integral operator of order Hs of \*on\$, % s ! is defined as:

**328** 

$$f\_s^a f(\mathbf{x}) = \begin{cases} \frac{1}{\Gamma(a)} \int\_s^\chi (\mathbf{x} - t)^{a-1} f(t) dt, & \mathbf{x} > t > s \ge 0, a > 0\\ f(\mathbf{x}), & a = 0 \end{cases} \tag{1}$$

**Definition 2.3:** The Caputo fractional derivative of order H of \*on : ;ok is defined as:

$$D\_s^\alpha f(\mathbf{x}) = \begin{cases} f\_s^{n-a} f^{(n)}(\mathbf{x}), & \mathbf{x} > \mathbf{s} \ge \mathbf{0}, n-1 < a < n\\ d^n f(\mathbf{x}), & a = n \end{cases} \tag{2}$$

On the one hand, for some certain properties of the operator "j F, it is obvious that when w ! 5 N s , and nol, we have "j F 5Nx  y x( y x(F 5NxF and "j Fn. On the other hand, properties of the operator tj F can be summarized shortly in the form of the following: for \*on\$ % s !, H s , nol, and w s !, we have tj Fn  <sup>z</sup> y F( 5NF, tj Ftj \* 5  tj F(\* 5  tj tj F\* 5, and tj F 5Nx  y x( y F(x( 5NF(x.

**Theorem 2.1:** If ;! H;, \*on\$ : , ;ok, and % s !, then "j Ftj F\* 5  \* 5 and tj F"j F\* 5  \* 5 { \*<sup>A</sup> <sup>N</sup>( vj<sup>|</sup> Ab : A01 , where 5Ns.

#### **3. General Form of Generalized Taylor's Formula**

In this section, we will introduce general form of generalized Taylor's formula that contains the Caputo definition for fractional derivatives. In fact, we need this generalization in the application of the RPS method in order to predict and find the multiplicity of solutions.

We will begin with the following definition which is needed throughout this work, especially, in the two succeeding sections. After that, we present a new and a fundamental theorem called general form of generalized Taylor's formula, which can formulate any function with certain properties in term of its fractional power series (FPS) representation.

**Definition 3.1:** A power series of the form

$$\sum\_{n=0}^{a} \sum\_{k=0}^{m-1} c\_{nk} \{t - t\_0\} \,^{k+na} \,, 0 \le m - 1 < a \le m \,, t \ge t\_0 \tag{3}$$

is called FPS about 1, where< is a variable and B:/'s are constants called the coefficients of the series.

As a special case, when <sup>1</sup>  , the expansion { { B:/ < J /(:F /01 . :01 is called a fractional Maclaurin series. Notice that in writing out the term corresponding to ; and 3 in Equation (3) we have adopted the convention that 1<sup>1</sup>  ! even when <sup>1</sup>. Also, when <sup>1</sup> each of terms of Equation (3) are vanished for ;}~3} and so. On the other hand, the FPS representation of Equation (3) always converges when < <sup>1</sup>. In the next lemma by " :F we mean that " <sup>F</sup> " <sup>F</sup> " <sup>F</sup> (;-times).

**Lemma 3.1:** Suppose that \* o n) 1 <sup>1</sup> & and " AF\* o n 1 <sup>1</sup> & for  ! 8 ; & ! where K! HK. Then we have:

$$\mathbf{329}$$

$$J\_{t\_0}^{(n+1)a} D\_{t\_0}^{(n+1)a} f(t) = \frac{\left(D\_{t\_0}^{(n+1)a} f\right)(\zeta)}{\Gamma\left( (n+1)a + 1 \right)} (t - t\_0)^{(n+1)a} \tag{4}$$

with <sup>1</sup> <sup>1</sup> & .

**Proof:** From the definition of the operator t F and by using the second mean value theorem for fractional integrals, one can find:

$$\begin{split} f\_{t\_0}^{(n+1)a} D\_{t\_0}^{(n+1)a} f(t) &= \frac{1}{\Gamma((n+1)a)} \int\_{t\_0}^t (t-\tau)^{(n+1)a-1} D\_{t\_0}^{(n+1)a} f(\tau) d\tau \\ &= \frac{\binom{D\_{t\_0}^{(n+1)a} f}{\Gamma((n+1)a)} \binom{\ell}{\ell}}{\Gamma((n+1)a)} \int\_{t\_0}^t (t-\tau)^{(n+1)a-1} d\tau \\ &= \frac{\binom{D\_{t\_0}^{(n+1)a} f}{\Gamma((n+1)a+1)} (t-t\_0)^{(n+1)a}}{\Gamma((n+1)a+1)} \end{split} \tag{5}$$

**Theorem 3.1:** Suppose that<\* o n) 1 <sup>1</sup> & , " AF\* o n 1 <sup>1</sup> & , and " AF\* can be differentiated K!-times on 1 <sup>1</sup> & for  !8 ; & !, where K! HK. Then:

$$f(t) = \sum\_{j=0}^{n} \sum\_{k=0}^{m-1} \frac{\binom{D^k D\_{t\_0}^{j\alpha} f}{\Gamma\{j\alpha + k + 1\}} (t - t\_0)^{k + j\alpha} + \frac{\binom{D\_{t\_0}^{(n+1)\alpha} f}{\Gamma\{(n+1)\alpha + 1\}} (t - t\_0)^{(n+1)\alpha}}{\Gamma\{(n+1)\alpha + 1\}} \tag{6}$$

with <sup>1</sup> <sup>1</sup> & .

**Proof:** From the certain properties of the operator t F , we have:

t :(F" :(F\*  t :F t <sup>F</sup> " <sup>F</sup> " :F\* t :F t J" J" :F\* <<<  t :F " :F\* - <sup>g</sup> / / " :F\*h 1 ( 3b J /01 1/  t :F" :F\* t :F 2- "/" :F\* 1 3b J /01 1/4  t :F t J" J" :F\* - "/" :F\* 1 m ;H & 3 & ! 1/(:F J /01 <<<<<<<<<  t :F " :F\* - "/" :F\* 1 3b J /01 1/ - "/" :F\* 1 m ;H & 3 & ! 1/(:F J< /01 (7)

On the other direction, if we keep the repeating of this process, then after ; -times of computations, we will obtain:

**330** 

$$f\_{\mathbf{t}\_0}^{\{n+1\}a} D\_{\mathbf{t}\_0}^{\{n+1\}a} f(\mathbf{t}) = f(\mathbf{t}) - \sum\_{j=0}^n \sum\_{k=0}^{m-1} \frac{\binom{D}{j} D\_{\mathbf{t}\_0}^{\{n\}} f(\mathbf{t}\_0)}{\Gamma(ja+k+1)} (\mathbf{t} - \mathbf{t}\_0)^{k+j\alpha}, \mathbf{t}\_0 \le \mathbf{t} \le \mathbf{t}\_0 + R \tag{8}$$

Thus, by using Lemma 3.1, the proof of the theorem will be complete.

**Remark 3.1:** We mention here that, if we fixed K!, then the series representation of \* in Theorem 3.1 will leads to the following expansion of \* about 1:

$$f(t) = \sum\_{j=0}^{n} \frac{D\_{t\_0}^{j\alpha} f(t\_0)}{\Gamma(ja+1)} (t - t\_0)^{j\alpha} + \frac{\left(D\_{t\_0}^{(n+1)\alpha} f\right)(\zeta)}{\Gamma\left((n+1)a+1\right)} (t - t\_0)^{(n+1)a} \tag{9}$$

with <sup>1</sup> <sup>1</sup> & , which is the same as of Generalized Taylor's series that obtained in [44] for H!.

As with any convergent series, this means that \* is the limit of the sequence of partial sums. In the case of general form of generalized Taylor's series, the partial sums are given as :  { { | y AF(/( J /01 1 : /(AF A01 . In general, \* is the sum of its general form of generalized Taylor's series if \*  d:. : . But on the other aspect as well, if we set :  \* : , then : is the remainder of general form of generalized Taylor's series. That is, : y :(F( 1 :(F <sup>1</sup> <sup>1</sup> & .

**Corollary 3.1:** If " :(F\* on <sup>1</sup> , where K HK!, then the reminder : of general form of generalized Taylor's series satisfies the inequality:

$$|R\_n(t)| \le \frac{M}{\Gamma((n+1)a+1)}(t-t\_0)^{(n+1)a}, t\_0 \le t \le d\tag{10}$$

#### **4. RPS Method for BVPs of Fractional Order**

In this section, we predict and find multiple solutions for BVPs of fractional order that admit multiple solutions by substituting a FPS expansion with undetermined coefficients through the given equation. From the FDE a recursion formula for the computation of the coefficients was derived. On the other hand, the coefficients in the FPS expansion can be computed recursively by differentiating the FDEs.

For convenience, the reader is referred to [28–30] in order to know more details about the classical RPS methods, including their construction, their motivation for use, their characteristics compared to the conventional method, and their applications for solving different categories of linear and nonlinear differential equations of different types and orders.

In fact, the main goal of our work is to predict and find out multiple series solutions for nonlinear BVPs of fractional order. To illustrate the basic idea of the RPS method for solving fractional BVPs analytically, we first consider the following nonlinear fractional functional equation:

$$D\_{r\_0}^{\alpha}u(r) = \mathcal{N}[u(r)], r \in \Omega, m - 1 < a \le m \tag{11}$$

subject to the boundary conditions:

$$\mathcal{B}\left(u(r), \frac{\partial u(r)}{\partial n}\right) = 0, r \in \Pi \tag{12}$$

**331** 

where " <sup>F</sup> is a Caputo fractional derivative of order ¤, is general nonlinear operator has Caputo fractional derivative term, is boundary operator, and £ is the boundary of the domain .

The crucial step of the RPS method for solving Equations (11) and (12) analytically is based on the fact that the boundary conditions of Equation (12) should be transcribed into equivalent form, so that, the new version conditions involves an unknown parameter so-called prescribed parameter ¥ and are split to:

$$\mathcal{B}^\*\left(u(r), \delta, \frac{\partial u(r)}{\partial n}\right) = 0, r \in \Pi, u(\alpha) = \beta \tag{13}$$

where @ H  is the forcing condition that comes from original conditions of Equation (12) and is the remainder boundary operator which contains the prescribed parameter ¥. Now, we investigate and apply the RPS method on the following problem:

$$D\_{r\_0}^{\alpha}u(r) = \mathcal{N}[u(r)], r \in \Omega, m - 1 < a \le m \tag{14}$$

subject to the split conditions:

$$\mathcal{B}^\* \left( u(r), \delta, \frac{\partial u(r)}{\partial n} \right) = 0, r \in \Pi \tag{15}$$

In order to apply the RPS method in the fractional sense, we assume firstly that we can apply the operator " =F" <sup>A</sup>  ! K ! ¦  ! on the term of )@ , in Equation (14) and also we suppose that all solutions @ that satisfy Equations (14) and (15) can be expanded by FPS representation as follows:

$$u(r) = \sum\_{l=0}^{\infty} \sum\_{j=0}^{m-1} c\_{lj} (r - r\_0)^{j + l\alpha}, r \ge r\_0 \tag{16}$$

where B=A are coefficients to be determined and 1 is the initial value of the independent variable of Equations (14) and (15).

Assume that Equation (14) satisfies the initial conditions @<sup>=</sup> 1  ¥= ¦  ! 8 K !, with ¥= as a prescribed parameters, where the unknown parameters ¥= can be determined later by substituting @ H  and/or any other constraint conditions into the obtained solution form. It is worth noting that some of ¥= may be known from the given initial conditions in Equation (15). As the first step in the prediction of multiple solutions, we set @11 be as an initial guess approximation of exact function solution @ of Equation (14). On the other hand, @11 will be of the form @11  @11 ' ¥1 ¥¥J which satisfies the known initial conditions in Equation (15) automatically. It was proved in [28–30] that the coefficients B1A in Equation (16) will take the form B1A  §| Ab  ! 8 K !. Therefore, we can consider the expansion formula:

$$u\_{0,0}(r) = \sum\_{l=0}^{m-1} \frac{\delta\_l}{j!} (r - r\_0)^l, r \ge r\_0 \tag{17}$$

as an initial guess approximation for solution of Equations (14) and (15). But on the other aspect as well, depending on Equations (16) and (17), we can write:

$$u(r) = \sum\_{j=0}^{m-1} \frac{\delta\_j}{j!} (r - r\_0)^j + \sum\_{l=1}^{\infty} \sum\_{j=0}^{m-1} c\_{lj} (r - r\_0)^{l+la}, r \ge r\_0 \tag{18}$$

Again, as the second step in the prediction of multiple solutions, we will let the double 3 ¨ to denote the 3 ¨-truncated series approximation of @ . That is:

$$u\_{k,l}(r) = \sum\_{j=0}^{m-1} \frac{\delta\_j}{j!} (r - r\_0)^j + \sum\_{l=1}^{k} \sum\_{j=0}^{l} c\_{lj} (r - r\_0)^{l+la}, r \ge r\_0 \tag{19}$$

where the indices counter ¨ and 3 whenever used mean that 3  ! 8 V and ¨  ! 8 K !.

Prior to applying the RPS technique for finding the values of coefficients B=A in the series expansion of Equation (19), we must define the residual function concept for the main nonlinear fractional functional Equation (14) as ª«¬  " <sup>F</sup> @ )@ , s1 and the following truncated 3 ¨-resudial function:

$$\operatorname{Res}\_{\{k,l\}}(r) = D\_{r\_0}^{\alpha} u\_{k,l}(r) - \mathcal{N}\left[u\_{k,l}(r)\right], r \ge r\_0 \tag{20}$$

As in [28–30], it is clear that ª«¬  for each o )1 1 & , where is the radius of convergence of Equation (18). In fact, this shows that " =F" <sup>A</sup> ª«¬  for each ¦  ! 8 V 3 and  ! 8 ¨, since the fractional derivative of a constant function in the Caputo sense is zero. In the mean time, the fractional derivatives " =F" <sup>A</sup> for each ¦  ! 8 V 3 and  ! 8 ¨ of ª«¬ and ª«¬ =A are matching at 1 ; it is obvious that:

$$D\_{r\_0}^{\{i-1\}\alpha} D\_{r\_0}^j \operatorname{Res} \{ r\_0 \} = D\_{r\_0}^{\{i-1\}\alpha} D\_{r\_0}^j \operatorname{Res}\_{\{l, \emptyset\}} \{ r\_0 \} = 0, l = 1, 2, 3, \dots, k, j = 0, 1, 2, \dots, l \qquad (21)$$

To obtain the value of coefficients B¯° in Equation (19) for ±  ! 8 V 3 and L  ! 8 ¨ , we apply the following subroutine: substitute ± L -truncated series approximation of @ into Equation (20), find the fractional derivative formula " ¯F" ° of ª«¬ ¯° at 1 , and then finally solve the obtained algebraic equation to get the required coefficients.

To summarize the computation process of RPS method in numerical values, we apply the following: fixed ¦! and run the counter  ! 8 ¨ to find ! -truncated series expansion of suggested solution, next fixed<¦  8 and run the counter  ! 8 ¨ to obtain the 8 -truncated series, and so on. In fact, to get ! -truncated series expansion for Equations (14) and (15), we use Equation (19) and write:

$$u\_{1,0}(r) = \delta\_0 + \delta\_1 \{r - r\_0\} + \dots + \delta\_{m-1} \{r - r\_0\}^{(m-1)} + c\_{10} \{r - r\_0\}^a \tag{22}$$

On the other hand, to determine the value of first unknown coefficient, B1, in Equation (22), we should substitute Equation (22) into both sides of the !-residual function that obtained from Equation (20), to get the following result:

$$\begin{aligned} \text{Res}\_{\{1,0\}}(r) &= c\_{1,0} \Gamma(a+1) \\ &- \mathcal{N} \left[ \delta\_0 + \delta\_1 (r - r\_0) + \dots + \delta\_{m-1} (r - r\_0)^{\{m-1\}} + c\_{10} (r - r\_0)^a \right] \end{aligned} \tag{23}$$

Now, depending on the result of Equation (21) for ¦  ! , Equation (23) gives B1  )§, y F( . Hence, using the !-truncated series expansion of Equation (22), the !-RPS approximation for Equations (14) and (15) can be expressed as:

$$u\_{1,0}(r) = \delta\_0 + \delta\_1(r - r\_0) + \dots + \delta\_{m-1}(r - r\_0)^{(m-1)} + \frac{\mathcal{N}[\delta\_0]}{\Gamma\{a+1\}}(r - r\_0)^a \tag{24}$$

Similarly, to find !!-truncated series expansion for Equations (14) and (15), we use Equation (19) and write:

$$\begin{split} u\_{1,1}(r) &= \delta\_0 + \delta\_1 (r - r\_0) + \dots + \delta\_{m-1} (r - r\_0)^{\{m-1\}} + \frac{N \lfloor \mathfrak{b}\_0 \rfloor}{\Gamma \{a+1\}} (r - r\_0)^a \\ &+ c\_{11} (r - r\_0)^{1+a} \end{split} \tag{25}$$

Again, to find out the value of second unknown coefficient, *c*11 in Equation (25), we must find and formulate !!-residual function based on Equation (20) and then substitute the form of @ in Equation (25) to find new discretized form of this residual function as follows:

$$\begin{aligned} \text{Res}\_{\{1,1\}}(r) &= D\_{r\_0}^a u\_{1,1}(r) - \mathcal{N}[u\_{1,1}(r)] \\ &= c\_{10} \Gamma(a+1) + c\_{11} \Gamma(a+2)(r - r\_0) \\ &- \mathcal{N} \left[ \delta\_0 + \delta\_1 (r - r\_0) + \dots + \delta\_{m-1} (r - r\_0)^{(m-1)} + \frac{\mathcal{N}[\delta\_0]}{\Gamma(a+1)} (r - r\_0)^a \right. \\ &+ c\_{11} (r - r\_0)^{1+a} \right] \end{aligned} \tag{26}$$

while, on the other hand, by considering Equation (20) for ¦  !! and applying the operator " to the both side of Equation (26), we get:

$$\begin{aligned} D\_{r\_0} \text{Res}\_{\{1,1\}}(r\_0) \\ &= c\_{1,1} \Gamma(a+2) \\ &- D\_{r\_0} \mathcal{N} \left[ \delta\_0 + \delta\_1 (r - r\_0) + \dots + \delta\_{m-1} (r - r\_0)^{\{m-1\}} \\ &+ \frac{\mathcal{N}[\delta\_0]}{\Gamma(a+1)} (r - r\_0)^a + c\_{11} (r - r\_0)^{1+a} \right]\_{r = r\_0} \end{aligned} \tag{27}$$

Now, using the fact that "ª«¬ 1  , we can easily obtain:

$$\begin{split} c\_{11} &= \frac{1}{\Gamma(a+2)} D\_{r\_0} \mathcal{N} \left[ \delta\_0 + \delta\_1 (r - r\_0) + \dots + \delta\_{m-1} (r - r\_0)^{(m-1)} \\ &+ \frac{\mathcal{N} \left[ \delta\_0 \right]}{\Gamma(a+1)} (r - r\_0)^a + c\_{11} (r - r\_0)^{1+a} \right]\_{r=r\_0} \end{split} \tag{28}$$

Hence, using the !! -truncated series expansion of Equation (25), the !! -RPS approximation for Equations (14) and (15) can be expressed as:

$$\begin{split} u\_{1,1}(r) &= \sum\_{j=0}^{m-1} \frac{\delta\_j}{j!} (r - r\_0)^j + \frac{N[\delta\_0]}{\Gamma(a+1)} (r - r\_0)^a \\ &\quad + \frac{1}{\Gamma(a+2)} D\_{r\_0} \mathcal{N} \left[ \delta\_0 + \delta\_1 (r - r\_0) + \dots + \delta\_{m-1} (r - r\_0)^{(m-1)} \right. \\ &\quad + \frac{\mathcal{N}[\delta\_0]}{\Gamma(a+1)} (r - r\_0)^a + c\_{11} (r - r\_0)^{1+a} \Big|\_{r = r\_0} (r - r\_0)^{1+a}, r \ge r\_0 \end{split} \tag{29}$$

This procedure can be repeated till the arbitrary order coefficients of FPS solution for Equations (14) and (15) are obtained. Moreover, higher accuracy can be achieved by evaluating more components.

**Remark 4.1:** It is worth indicating that there are still unknown prescribed parameters ¥= in the series expansion of Equation (19) (and simply in Equation (29)) that should be determined. It is essential that the existence of a unique or multiple solutions in terms of Equation (19) (and simply in Equation (29)) for the original BVP which is covered by Equations (11) and (12) depends on the fact that whether the forcing condition @ H<  and/or any other constraint condition in Equation (13) admits a unique or multiple values for the formally introduced prescribed parameters ¥=. This stage is called rule of multiplicity of solutions that is a criterion in order to know how many solutions the BVP in Equations (11) and (12) admits.

Anyhow, as the final step in the construction, if we substitute @ H  and/or any other constraint conditions into the obtained solution form of Equation (19) (and simply in Equation (29)), then we obtain a system of nonlinear algebraic equations in the prescribed variables ¥1 ¥ ¥7¥J (here, we must recall that some of ¥= ¦  ! K ! may be known from Equation (15)) which can be easy solved using symbolic computation software such as MAPLE !V or MATHEMATICA ZC. In fact, if we substitute these values of prescribed parameters in the obtained solution form of Equation (19) (and simply in Equation (29)), then discretized form of the 3 ¨-truncated series approximation of @ of Equations (11) and (12) (and simply !!-truncated series approximation of @ as given by Equation (26)) will be obtained.

#### **5. Applications and Numerical Discussions**

The application problems are carried out using the proposed RPS method, which is one of the modern analytical techniques because of its iterative nature; it can handle any kind of boundary conditions and other constraints. The RPS method doesn't have mathematical requirements about the multiple solutions of fractional BVPs to be solved; the RPS method is also very effective in identifying global predicted solutions, and provides a great flexibility in choosing the initial guess approximations. However, in order to verify the computational efficiency of the designed RPS method, two nonlinear models are performed, one of them arises in mixed convection flows and the other one arises in heat transfer, which both admit multiple solutions. In the process of computation, all the symbolic and numerical computations were performed by using the MATHEMATICA ZC software package.

Throughout this section, we will try to give the results of the two applications; however, in some cases we will switch between the results obtained for the applications in order not to increase the length of the paper without the loss of generality for the remaining application and results. However, by easy calculations we can collect further results and discussion for the desire application.

**Application 5.1:** The aim of this application is to apply the RPS method to analyze a kind of model in mixed convection flows namely, combined forced and free flow in the fully developed region of a vertical channel with isothermal walls kept at the same temperature. In this model, the fluid properties are assumed to be constant and the viscous dissipation effect is taken into account. The set of governing balance equations for the velocity field is reduced to the following [45,46]:

$$D\_0^\alpha u(\mathbf{y}) = \frac{\Psi}{16} \left(\frac{du(\mathbf{y})}{d\mathbf{y}}\right)^2, 3 < a \le 4, 0 \le \mathbf{y} \le 1\tag{30}$$

subject to boundary conditions:

$$u'(0) = u^{\prime\prime\prime}(0) = u(1) = 0, \int\_0^1 u(\mathbf{y}) \, d\mathbf{y} = 1 \tag{31}$$

where @ and I are dimensionless velocity and transversal coordinate, respectively, and also @  · ·¸ , I  ¹ <sup>º</sup> , »O  Rº¼ ½¾ , ^  \$½¾ / , O  Rº·¸ ¯ , and µ  ¿«ÀÁª« in which ÂJ Ã Ä + Bq % 3 ± ¿« ÀÁ, and ª« are mean 5-component of the fluid velocity, fluid velocity, channel half-width, acceleration due to gravity, coefficient of thermal expansion, specific heat at constant pressure, dynamic viscosity, thermal conductivity, kinematic viscosity, Gebhart number, Prandtl number, and Reynolds number, respectively.

Next, we will show how one can find out the existence of multiple solutions for Equations (30) and (31) in aforesaid range for < I !. To do so, we consider firstly Equations (30) and (31) and suppose that @  ¥1 and @¶¶  ¥7. So, Equations (30) and (31) can be modified into the following form:

$$D\_0^\alpha u(\mathbf{y}) = \frac{\Psi}{16} \left(\frac{du(\mathbf{y})}{d\mathbf{y}}\right)^2, 3 < a \le 4, 0 \le \mathbf{y} \le 1\tag{32}$$

subject to the split conditions:

$$u(0) = \delta\_0, u'(0) = 0, u''(0) = \delta\_2, u''''(0) = 0\tag{33}$$

where @ !  is the additional forcing condition and <sup>Å</sup> @ I <sup>1</sup> I  ! is the additional constraint condition.

Now, we apply the RPS method on Equations (32) and (33), where prescribed parameters ¥1 and ¥7 which are played an important and fundamental role to realize about multiplicity of solutions, will be obtained later by substituting the additional forcing condition and the additional constraint condition in resulting expansion formula that approximate Equations (32) and (33).

According to Equation (18), we assume that the series solution of Equations (32) and (33) can be written as:

$$u(\mathbf{y}) = \sum\_{j=0}^{3} \frac{\delta\_j}{j!} \mathbf{y}^j + \sum\_{l=1}^{\infty} \sum\_{j=0}^{3} c\_{lj} \mathbf{y}^{j+la} \tag{34}$$

where ¥  @¶  and ¥M  @¶¶¶  which are hold from the conditions of Equation (33). Therefore, the initial guess approximation can be constructing as @11 I  ¥1 & ÆÇ <sup>7</sup> <sup>I</sup><sup>7</sup> . Next, according to Equations (19) and (20) the 3 ¨-truncated series approximation of @ I and the 3 ¨-residual function of Equation (32) can be defined and thus constructed, respectively, as:

$$\operatorname{Res}\_{\{k,l\}}\{\mathbf{y}\} = D\_0^\alpha u\_{k,l}(\mathbf{y}) - \frac{\Psi}{16} \left(\frac{du\_{k,l}(\mathbf{y})}{dy}\right)^2 \tag{35}$$

$$u\_{k,l}(\mathbf{y}) = \delta\_0 + \frac{\delta\_2}{2}\mathbf{y}^2 + \sum\_{l=1}^{k} \sum\_{j=0}^{l} c\_{lj}\mathbf{y}^{j+la} \tag{36}$$

However, to determine the value of coefficient B1 , we find out (1,0)-truncated series approximation of @ I as @1 I  ¥1 & §Ç <sup>7</sup> <sup>I</sup><sup>7</sup> & B1I<sup>F</sup> and (1,0)-residual function as ª«¬ 1 I  B1m H&! <sup>È</sup> [ ¥7I & HB1IF7 . On the other aspect as well, by using Equation (21) for ¦  ! and substituting I, we obtain B1  .

Similarly, to find the value of coefficient B, we evaluate !!-truncated series approximation of @ I as @ I  ¥1 & ÆÇ <sup>7</sup> <sup>I</sup><sup>7</sup> & BI(F and !! -residual function as ª«¬ I  Bm H&8I <sup>È</sup> [ ¥7I & !&HBIF7 . Thus, for ¦  !! , we conclude that "1ª«¬ I  Bm H&8 <sup>È</sup> <sup>T</sup> ¥7 & H !&HBIF ¥7I & !&HBIF , while the substitution of I leads to B  .

To evaluate the value of coefficient B7, we need to write @7 I  ¥1 & ÆÇ <sup>7</sup> <sup>I</sup><sup>7</sup> & B7I7(F and ª«¬ 7 I  B7 y F(M <sup>7</sup> <sup>I</sup><sup>7</sup> <sup>È</sup> [ ¥7I&B7 8&HI(F7. However, by considering the fact that "1 7ª«¬ 7  , we can easily find B7  ÈÆÇ Ç Ty M(F . Similarly, the continuation in the same manner will leads also to BM  . According to the initial guess approximation and the form of terms in Equation (32) taking into account the form of Equations (35) and (36), we can conclude that B=A  for  !V. Therefore, according to Equation (34) the FPS solution of Equations (32) and (33) can be written in the form of the following expansion:

$$u(\mathbf{y}) = \delta\_0 + \frac{\delta\_2}{2}\mathbf{y}^2 + \frac{\Psi \delta\_2^2}{8\Gamma(3+a)}\mathbf{y}^{2+a} + \sum\_{l=2}^{\infty} c\_{l2}\mathbf{y}^{2+la} \tag{37}$$

and hence the 3 8-truncated series approximation of @ I can reformulated as:

$$u\_{k,2}(\mathbf{y}) = \delta\_0 + \frac{\delta\_2}{2}\mathbf{y}^2 + \frac{\Psi \delta\_2^2}{8\Gamma(3+\alpha)}\mathbf{y}^{2+\alpha} + \sum\_{l=2}^{k} c\_{l2}\mathbf{y}^{2+l\alpha} \tag{38}$$

Again, to determine the value of coefficient B77, we need solve the equation "1 F"1 7ª«¬ 77  which gives B77  7(FÈÇÆÇ É [Ry M(7F . Similarly, we have BM7  8&HµM¥7 R Pm V&H7 & PHm V & H8&8mV&8H&HmV&8HQ!8PmV&H8mV&VH, and so on. Consequently, based on the structure of Equation (38) the V8-truncated series approximation of @ I generated from the RPS method can be written as:

$$\begin{aligned} &u\_{3,2}(\mathbf{y})\\ &=\delta\_0 + \frac{\delta\_2}{2}y^2 + \frac{\Psi \delta\_2^2}{8\Gamma(3+a)}y^{2+a} + \frac{(2+a)\Psi^2 \delta\_2^3}{64\Gamma(3+2a)}y^{2+2a} \\ &+\frac{(2+a)\Psi^3 \delta\_2^4 \{4\Gamma(3+a)^2 + 4a\Gamma(3+a)^2 + 2\Gamma(3+2a) + a\Gamma(3+2a)\}}{1024\Gamma(3+a)^2\Gamma(3+3a)}y^{2+3a} \end{aligned} \tag{39}$$

It is clear that, Equation (39) contain two unknown certain parameters which are ¥1 and ¥7. To determine their introductory values substituting the forcing condition @ !  and the constraint condition <sup>Å</sup> @ I <sup>1</sup> I  !, and finally selecting some numerical values of µ.

Now, to be specific, we consider two cases according to Equation (39) which consist of  8 and  8. On the other hand, we generate and obtain the !Z8-truncated series approximation of @ I using the same procedure. However, various values of ¥1 and ¥7 have been calculated and listed in Tables 1 and 2 when H  VCX and when HP, respectively. For simplicity and not to conflict, we will let @/© I to denote the first approximate solution of @ I and @/© <sup>7</sup> I to denote the second approximate solution of @ I.

Consequently, we conclude that the RPS method furnishes multiple solutions for Equations (30) and (31). It is worth mentioning here that when H  VCX, Table 1 indicates the existence of two solutions at  8, so that, @  !CP\X\WVPYZP, @¶¶  8CX\!8YPZ\PP for the first branch solution and @  !VCZ!W\!!8P8V, @¶¶  !P\CVYWPPZP8XW for the second branch solution. In fact, these results answer the question how many solutions the nonlinear BVP in Equations (30) and (31) admits? The same procedure has been done at the case  8. As we see from Table 1, there exist multiple solutions namely @  !CW!!8WWW!\8 @¶¶  VC!!XZY!PZV for the first branch solution and @  !WC8X8YYZX\YW @¶¶  !VXCZYPYP8W!VW for the second branch solution. Similar conclusion can be achieved when HP as shown in Table 2.

**Table 1.** The approximate numerical values of ¥1 / and ¥7 /, 3  !8 at  8 and  8 when H  VCX for @S7 I.


**Table 2.** The approximate numerical values of ¥1 / and ¥7 /, 3  !8 at  8 and  8 when HP for @S7 I.


The RPS technique has an advantage that it is possible to pick any point in the interval of integration and as well the approximate multiple solutions and all their fractional derivatives will be applicable. In other words, continuous approximate solutions can be obtained. Our next goal, is to show the mathematical behavior of the obtained multiple solutions geometrically. To do so, we plot the first and the second solutions obtained from the !Z8-truncated series approximation of @ I at µ  8 and µ  8 when<H  VCX in Figure 1, while in Figure 2 we depict the first and the second approximate solutions at the same values when HP.

The effective calculations of the two approximate branches solutions for Equations (30) and (31) with respect to some certain specific values of µ on @ and @¶¶ is explored next in which the obtained results are generated from the !Z8-truncated series approximation of @ I. Table 3 gives the effect of the numeric value of µ when H  VCX, while Table 4 gives the effect of the numeric value of µ when HP. The numeric value of µ lie within the range )\\, in step of 8. It is to be noted that, when the values of µ increasing gradually within mentioned range, the value of @¶¶ decreasing as well as the value of @ increasing for both branch solutions and for both order of derivatives.

We mention here that, the case of µ correspond either to a very small viscous dissipation heating or to negligible buoyancy effects. However, Equations (30) and (31) are easily solved and admit the unique solution @ I  <sup>M</sup> <sup>7</sup> !I7 for both H  VCX and HP . On the other direction, from the last tables, it can be seen that our results of the RPS method agree best with method of [34] when H  VCX and method of [35] when HP.

**Figure 1.** Multiple solutions of Equations (30) and (31) when H  VCX: @S7 I: red color, @S7 <sup>7</sup> I: blue color at (**a**) µ  8 and (**b**) µ  8.

**Figure 2.** Multiple solutions of Equations (30) and (31) when HP: @S7 I: red color, @S7 <sup>7</sup> I: blue color at (**a**) µ  8 and (**b**) µ  8.


**Table 3.** The effect values of µ on ¥1 / and ¥7 /, 3  ! 8 for the first and the second branch solutions when H  VCX.

**Table 4.** The effect values of µ on ¥1 / and ¥7 /, 3  !8 for the first and the second branch solutions when HP.


To measure the accuracy and the efficiency of the proposed RPS method for predicting and finding the multiple solutions for Equations (30) and (31), we report the residual error function at  8 and  8 when H  VCX and when HP in Tables 5 and 6, respectively, in which the obtained results are generated from the !Z8-truncated series approximation of @ I. The residual error function is defined using Equation (35) in which the grid points are building as I=  1 ¦, ¦  ! 8 !. For simplicity and not to conflict, we will let ª«¬ /© I to denote the residual error function of the first approximate solution @/© I of @ I and ª«¬ /© <sup>7</sup> I to denote the residual error function of the second approximate solution @/© <sup>7</sup> I of @ I.

In fact, the residual errors measure the extent of agreement between the !Z8 th-order approximate RPS solutions and unknowns closed form solutions which are inapplicable in general for such nonlinear equations. However, from the tables, it can be seen that the RPS technique provides us with the accurate approximate solutions and explicates the rapid convergence in approximating the multiple solutions for Equations (30) and (31).

**Table 5.** The values of absolute residual error function ª«¬ S7 <sup>Ñ</sup> I, 3  !8 at  8 and  8 of @S7 I when H  VCX.


**Table 6.** The values of absolute residual error function ª«¬ S7 <sup>Ñ</sup> I, 3  !8 at  8 and  8 of @S7 I when HP.


In Tables 7 and 8 we tabulate the values of the approximate multiple solutions at the final grid node I! that generated from the !Z8-truncated series approximation of @ I at  8 and  8 when H  VCX and when HP. In fact, we do this to facilitate the calculations in order to show the validity and accuracy of the proposed RPS method in predicting and finding the multiple approximate solutions. In these tables, we can find that the values of @S7 / !, 3  ! 8 agree nicely and completely the forcing condition @ !  and the constraint condition Å @S7 / I <sup>1</sup> , 3  ! 8.

**Table 7.** The approximate value of forcing condition @S7 / ! and constraint condition Å @S7 / I <sup>1</sup> I, 3  !8 at  8 and  8 when H  VCX.


**Table 8.** The approximate value of forcing condition @S7 / ! and constraint condition Å @S7 / I <sup>1</sup> I, 3  !8 at  8 and  8 when HP.


In order to study the behavior of multiple approximate solutions in a better view, we plot the normalized of the two branch approximate solutions of Equations (30) and (31) with respect to @S7 / , 3  !8 at some specific values of and H in which the obtained results are generated from the !Z8-truncated series approximation of @ I. However, Figure 3 shows the normalized function ÜÝÇ Þ ÜÝÇ 1 , 3  !8 at  8 and  8 when H  VCX, while Figure 4 shows the normalized function ÜÝÇ Þ ÜÝÇ 1 , 3  !8 at  8 and  8 when HP. In these figures, we can see the almost similarity in the behavior of the two branches approximate solutions at the two mentioned specific value of when H  VCX and when HP.

**Figure 3.** Multiple solutions of Equations (30) and (31) via dimensionless transversal coordinate I when H  VCX: ÜÝÇ Þ ÜÝÇ 1 : red color, ÜÝÇ <sup>Ç</sup> Þ ÜÝÇ <sup>Ç</sup> 1 : blue color at (**a**) µ  8 and (**b**) µ  8.

**Figure 4.** Multiple solutions of Equations (30) and (31) via dimensionless transversal coordinate I when HP: ÜÝÇ Þ ÜÝÇ 1 : red color, ÜÝÇ <sup>Ç</sup> Þ ÜÝÇ <sup>Ç</sup> 1 : blue color at (**a**) µ  8 and (**b**) µ  8.

**Application 5.2:** Fins are extensively used to enhance the heat transfer between a solid surface and its convective, radiative, or convective radiative surface. Finned surfaces are widely used, for instance, for cooling electric transformers, the cylinders of aircraft engines, and other heat transfer equipment. The temperature distribution of a straight rectangular fin with a power-law temperature dependent surface heat flux can be determined by the solutions of a one-dimensional steady state heat conduction equation which, in dimensionless form, is given as follows [47,48]:

$$D\_0^{\alpha} \theta(\varkappa) \theta^3(\varkappa) = \frac{4}{25}, 1 < \alpha \le 2, 0 \le \varkappa \le 1\tag{40}$$

subject to the boundary conditions:

$$
\theta'(0) = 0, \theta(1) = 1 \tag{41}
$$

The prediction and construction of multiple solutions for BVPs of fractional order is the fundamental target of this paper. Next, we will show in brief steps and calculations how we can predict and find out existence of multiple solutions for Equations (40) and (41). To do so, we consider firstly, Equations (40) and (41) and suppose that ß  ¥1. So, a new discretized form of Equations (40) and (41) can be obtained as follows:

$$D\_0^{\alpha} \theta(\mathbf{x}) \theta^3(\mathbf{x}) = \frac{4}{25}, 1 < \alpha \le 2, 0 \le \mathbf{x} \le 1 \tag{42}$$

subject to the split conditions:

$$
\theta(0) = \delta\_0, \theta'(0) = 0 \tag{43}
$$

where ß !  ! is the additional forcing condition. Here, ¥1 denotes temperature of the fin tip and will be determined later by the rule of multiplicity of solutions from the process of computations thought the RPS technique.

Similar to the previous procedure and discussions that used in Application 5.1, the FPS solution and the residual function for Equations (42) and (43) will take, respectively, the following form:

$$\theta(\mathbf{x}) = \delta\_0 + \sum\_{l=1}^{\infty} \sum\_{j=0}^{1} c\_{lj} \mathbf{x}^{j+l\alpha} \tag{44}$$

$$\text{Res}\{\mathbf{x}\} = D\_0^{\alpha} \theta(\mathbf{x}) \theta^3(\mathbf{x}) - \frac{4}{25} \tag{45}$$

while the 3 ¨-truncated series approximation of ß 5 and the 3 ¨-resudial function that are derived from Equations (44) and (45) can be formulated, respectively, in form of:

$$\theta\_{k,l}(\mathbf{x}) = \delta\_0 + \sum\_{l=1}^{K} \sum\_{j=0}^{l} c\_{lj} \mathbf{x}^{j+la} \tag{46}$$

$$\operatorname{Res}\_{\{k,l\}}\{\boldsymbol{x}\} = D\_0^{\alpha} \theta\_{k,l}(\boldsymbol{x}) \theta\_{k,l}^3(\boldsymbol{x}) - \frac{4}{25} \tag{47}$$

It is to be noted that the !-truncated series solution of Equations (42) and (43) is ß1 5  ¥1 & B15<sup>F</sup> and the ! -residual function is ª«¬ 1 5  B1m H&! ¥1 & B15FM <sup>R</sup> 7U . Thus, using Equation (21) for ¦  !, we get B1  <sup>R</sup> 7UÆ Éy F( . Similarly, the !!-truncated series solution is ß 5  ¥1 & <sup>R</sup> 7UÆ Éy F( 5<sup>F</sup> & B5(F and the !!-residual function is:

ª«¬ 5  !8B 8W¥1 & !8Bm 8&H 8W¥1m !&H ¢ 5(F & !8B 7 8W¥1 <sup>7</sup> & 8PB <sup>7</sup> m 8&H 8W¥1 7m !&H ¢ 57(7F <<<<<<<<<<<<<<<<<<<<<<<<<<<& PB M 8W¥1 <sup>M</sup> & !8B <sup>M</sup> m 8&H 8W¥1 Mm !&H ¢ 5M(MF & 8WY VXY8W¥1 7m !&HM <sup>5</sup>MF & !X8 !WY8W¥1 Tm !&H7 <sup>5</sup>7F & !X8B !WY8W¥1 Ûm !&H7 & YPBm 8&H !WY8W¥1 Ûm !&HM¢ 5(MF & P\ Y8W¥1 Rm !&H 5F & XYB Y8W¥1 Um !&H & P\Bm 8&H Y8W¥1 Um !&H7¢ 5(7F & P\B 7 Y8W¥1 [m !&H & P\B <sup>7</sup> m 8&H Y8W¥1 [m !&H7¢ 57(MF & ¥1 MBm 8&H5 & V¥1 7B <sup>7</sup> m 8&H57(F & V¥1B <sup>M</sup> m 8&H5M(7F & B <sup>R</sup> m 8&H5R(MF (48)

More precisely, according to Equation (21) the solution of equation "1ª«¬  will gives B  . Thus, based on the initial guess approximation and the form of terms of Equation (42) taking into account the form of Equations (46) and (47), it easy to see that B=  , ¦  ! 8 V . Therefore, according to Equation (44) a new discretized form of FPS solution for Equations (42) and (43) can be obtained and expressed as:

$$\theta(\mathbf{x}) = \delta\_0 + \sum\_{l=1}^{\infty} c\_{l0} \mathbf{x}^{l\alpha} \tag{49}$$

and hence the 3 -truncated series approximation of ß 5 can formulated as:

$$\theta\_{k,0}(\mathbf{x}) = \delta\_0 + \frac{4}{25\delta\_0^3 \Gamma(\alpha + 1)} \mathbf{x}^a + \sum\_{l=2}^{k} c\_{l0} \mathbf{x}^{la} \tag{50}$$

In the shape of shapes by continuing in this procedure and using Equations (46) and (47) taking into account Equation (21), we can easily obtained that B71  RT [7UÆÝy (7F , BM1  Û7 My (FÇ(7y (7F U[7UÆ y (FÇy (MF , and so on. Consequently, based on Equation (50) the V-truncated series approximation of ß 5 generated from the RPS method can be written as:

$$\begin{split} \theta\_{3,0}(\mathbf{x}) &= \delta\_0 + \frac{4}{25\delta\_0^3 \Gamma(1+a)} \mathbf{x}^a - \frac{48}{625\delta\_0^7 \Gamma(1+2a)} \mathbf{x}^{2a} \\ &+ \frac{192\{3\Gamma(1+a)^2 + 2\Gamma(1+2a)\}}{15625\delta\_0^{11} \Gamma(1+a)^2 \Gamma(1+3a)} \mathbf{x}^{3a} \end{split} \tag{51}$$

It is clear that, all the terms in Equation (51) contain an unknown certain parameter ¥1 and to determine its introductory values we must substitute the boundary condition ß ! 1 back into Equation (51) to obtain a nonlinear algebraic equation in one variable, which can be easy solved using symbolic computation software. But on the other aspect as well, if we generate and obtain

 V-truncated series approximation of ß 5 by using the same procedure discussed, then two various values of ¥1 have been calculated and listed in Table 9 when H  !CX and when H8.

**Table 9.** The approximate numerical values of ¥1 /, 3  !8 when H  !CX and when H8 for ßM111 5.


It is clear from the table that two ¥-plateaus can be identified and consequently we conclude that the RPS method furnishes multiple solutions to Equations (40) and (41). It is worth mentioning here that Equation (46) (and simply Equation (51)) indicates existence of two solutions. On the other hand, the existence of a unique or multiple solutions in terms of Equation (46) (and simply Equation (51)) for the original BVP which is covered by Equations (40) and (41) depends on the fact that whether the forcing condition ß !  ! admits a unique or multiple values for the formally introduced prescribed parameters ¥1.

Finally, in Figure 5 we plot the first and the second approximate multiple solutions of Equations (40) and (41) that obtained from the V-truncated series approximation of ß 5 when<H  !CX and when H8. In fact, we do this for the same reasons that mentioned in the Application 5.1, where the same conclusion can be obtained too. On the other direction, as in the previous application, we can see that the sketch of the two branches approximate RPS solutions that the problem admitted are agree best and nicely with method of [34] when H  !CX and method of [35] H8.

**Figure 5.** Multiple solutions of Equations (30) and (31): ßM111 5 : red color, ßM111 <sup>7</sup> 5: blue color when (**a**) H  !CX and (**b**) H8.

#### **6. Conclusions**

It is very important not to lose any solution of nonlinear FDEs with boundary conditions in engineering and physical sciences. In this regard, the present paper has introduced a new methodology namely the RPS method to prevent this, so that the presented method is not only able to predict the existence of multiple solutions, but also to calculate all branches of solutions effectively at the same time by using an appropriate initial guess approximation. We also noted that the RPS solutions were computed via a simple algorithm without any need for perturbation techniques, special transformations, or discretization. The validity of this method has been checked

by two nonlinear models, one of them arises in mixed convection flows and the other one arises in heat transfer, which both admit multiple or dual solutions.

### **Acknowledgments**

The authors would like to express their thanks to unknown referees for their careful reading and helpful comments.

#### **Conflicts of Interest**

The authors declare no conflict of interest.

#### **References**



Reprinted from *Entropy*. Cite as: El-Ajou, A.; Arqub, O.A.; Zhour, Z.A.; Momani, S. New Results on Fractional Power Series: Theories and Applications. *Entropy* **2013**, *15*, 5305-5323.

### *Article*
