Approximate Solutions of an Extended Multi-Order Boundary Value Problem by Implementing Two Numerical Algorithms

Abstract

In this paper, we establish several necessary conditions to confirm the uniqueness-existence of solutions to an extended multi-order finite-term fractional differential equation with double-order integral boundary conditions with respect to asymmetric operators by relying on the Banach’s fixed-point criterion. We validate our study by implementing two numerical schemes to handle some Riemann–Liouville fractional boundary value problems and obtain approximate series solutions that converge to the exact ones. In particular, we present several examples that illustrate the closeness of the approximate solutions to the exact solutions.

1. Introduction

Fractional calculus is extending quickly, and its interesting and attractive applications are perfectly utilized in different parts of science [1,2,3]. It has appeared in financial models [4], optimal control [5,6], chaotic systems [7], epidemiological models [8,9], engineering [10,11], etc. Particularly, the fractional systems of boundary value problems (FBVP) of fractional differential equations usually yield other operational mathematical models for the description of special chemical, physical, and biological processes, which one can find in recently published works [12,13,14,15,16,17,18,19]. Along with these real models describing the phenomena, many mathematicians conduct research on the existence theory of solutions for different abstract structures of FBVPs with general boundary conditions including three-point, multi-point, multi-order, multi-strip, and nonlocal integral ones [20,21,22,23,24,25,26,27,28,29].

Several studies have also concentrated on the numerical techniques to obtain the analytical and approximate solutions of FBVPs. New numerical methods are introduced by researchers that have improved the convergence rate and error resulting from the approximate solutions. Examples of these methods and how to use them are Haar wavelet method [30,31], CAS wavelet method [32], homotopy analysis transform method (HATM) [33], q-HATM [34], Bernstein polynomials [35], iterative reproducing kernel Hilbert space method [36], Legendre functions with fractional orders [37], variational iteration method [38], and so on.

Since multi-term multi-order fractional differential equations have appeared in a wide range of fields, many mathematicians have started to review the properties and numerical solutions of this type of fractional differential equations. On the other side, because most of the time the exact solution cannot be found or it is very difficult to find, various numerical techniques have been applied for such FBVPs to obtain the approximate solutions. For instance, Bolandtalat, Babolian and Jafari [39] compared the convergence effects of exact and numerical solutions of multi-order fractional differential equations by means of Boubaker polynomials. In 2016, Hesameddini, Rahimi, and Asadollahifard [40] presented a new version of the reliable algorithm to solve multi-order fractional differential equations and investigated the convergence of it. Firoozjaee et al. [41] implemented a numerical approach on a multi-order fractional differential equation with mixed boundary-initial conditions. Recently, Dabiri and Butcher [42] invoked a numerical technique based on the spectral collocation methods and obtained the numerical solutions of multi-order fractional differential equations subject to multiple delays.

In recent years, many FBVPs with integral boundary conditions have been formulated by researchers of this field. Ali, Sarwar, Zada, and Shah [43] developed some conditions with the aid of topological degree results for confirming the existence of solutions to the nonlinear integral FBVP

$$\left\{\begin{array}{c}{\displaystyle {}^c{\mathfrak{D}}_{0^{+}}^{\varrho }\upsilon \left(z\right)=\mathit{h}\left(z,\upsilon \left(\mu z\right)\right),z\in I:=[0,1],\mu \in (0,1),}\hfill \\ c_1\upsilon \left(0\right)+c_2\upsilon \left(1\right)={\mathfrak{I}}_{0^{+}}^{\varrho }{\phi }_1(1,\upsilon \left(1\right)),c_3{\upsilon }^{\prime }\left(0\right)+c_4{\upsilon }^{\prime }\left(1\right)={\mathfrak{I}}_{0^{+}}^{\varrho }{\phi }_2(1,\upsilon \left(1\right)),\hfill \end{array}\right.$$

in which $\varrho \in (1,2]$, $c_1,c_2,c_3,c_4\in {\mathbb{R}}^{+}$ and $h,{\phi }_1,{\phi }_2\in C(I\times \mathbb{R},\mathbb{R})$. ${}^c{\mathfrak{D}}_{0^{+}}^{\varrho }$ denotes the Caputo fractional derivative of order $\varrho$ and ${\mathfrak{I}}_{0^{+}}^{\varrho }$ is the Riemann–Liouville fractional integral of order $\varrho$. Liu, Li, Dai, and Liu [44] implemented the fixed point techniques to establish the existence and uniqueness of solutions for the nonlocal integral FBVP

$$\left\{\begin{array}{c}{\displaystyle {\mathfrak{D}}_{0^{+}}^{\varrho }\upsilon \left(z\right)+\psi \left(z\right)\mathit{h}\left(z,\upsilon \left(z\right)\right)=0,z\in (0,1),}\hfill \\ \upsilon \left(0\right)={\upsilon }^{\prime }\left(0\right)=\dots ={\upsilon }^{(k-2)}\left(0\right)=0,{\upsilon }^{\prime }\left(1\right)=p{\mathfrak{I}}_{0^{+}}^{\mu }\upsilon \left(\xi \right),\hfill \end{array}\right.$$

where $\varrho \in (k-1,k]$, $\xi \in (0,1]$, $p,\mu >0$, ${\displaystyle \frac{p\Gamma \left(\varrho \right){\xi }^{\varrho +\mu -1}}{\Gamma (\varrho +\mu )}}<1$ and ${\mathfrak{D}}_{0^{+}}^{\varrho }$ is the Riemann–Liouville fractional derivative of order $\varrho$. In 2018, Padhi, Graef, and Pati [45] studied positive solutions for the given fractional differential equation with Riemann–Stieltjes integral conditions

$$\left\{\begin{array}{c}{\displaystyle {\mathfrak{D}}_{0^{+}}^{\varrho }\upsilon \left(z\right)+\psi \left(z\right)\mathit{h}\left(z,\upsilon \left(z\right)\right)=0,z\in (0,1),}\hfill \\ {\displaystyle \upsilon \left(0\right)={\upsilon }^{\prime }\left(0\right)=\dots ={\upsilon }^{(k-2)}\left(0\right)=0,{\mathfrak{D}}_{0^{+}}^{\varpi }\upsilon \left(1\right)={\int }_0^1\phi (r,\upsilon \left(r\right))\mathrm{d}A\left(r\right),}\hfill \end{array}\right.$$

where $\varrho \in (k-1,k]$ with $k>2$ and $1\le \varpi \le \varrho -1$.

In 2021, Thabet, Etemad, and Rezapour [46] designed and discussed the notion of the existence for possible solutions of a coupled system of the Caputo conformable FBVPs of the pantograph differential equation by

$$\left\{\begin{array}{c}{}^{\mathcal{CC}}{\mathfrak{D}}_{z_0}^{\varrho ,{\sigma }_1^{*}}v\left(z\right)={\tilde{\mathcal{P}}}_1(z,m\left(z\right),m(\ell z)),z\in [z_0,\tilde{K}],z_0\ge 0,\hfill \\ {}^{\mathcal{CC}}{\mathfrak{D}}_{z_0}^{\varrho ,{\sigma }_2^{*}}m\left(z\right)={\tilde{\mathcal{P}}}_2(z,v\left(z\right),v(\ell z)),\hfill \end{array}\right.$$

with three-point RL-conformable integral conditions

$$\left\{\begin{array}{c}v\left(z_0\right)=0,c_{1}v\left(\tilde{K}\right)+c_2{}^{\mathcal{RC}}{\mathfrak{I}}_{z_0}^{\varrho ,{\theta }^{*}}v\left(\delta \right)=w_1^{*},\hfill \\ m\left(z_0\right)=0,c_1^{*}m\left(\tilde{K}\right)+c_2^{*}{}^{\mathcal{RC}}{\mathfrak{I}}_{z_0}^{\varrho ,{\theta }^{*}}m\left(\nu \right)=w_2^{*},\hfill \end{array}\right.$$

in which $\varrho \in (0,1]$, ${\sigma }_1^{*},{\sigma }_2^{*}\in (1,2)$, $\delta ,\nu \in (z_0,\tilde{K})$, $c_1,c_2,c_1^{*},c_2^{*},w_1^{*},w_2^{*}\in \mathbb{R}$, $\ell \in (0,1)$ and ${\tilde{\mathcal{P}}}_1,{\tilde{\mathcal{P}}}_2\in C([z_0,\tilde{K}]\times \mathbb{R}\times \mathbb{R},\mathbb{R})$. In all the above fractional models with integral conditions, only the required conditions of the existence of solutions have been investigated and FBVPs have not been solved numerically. Due to the complexity of the structure of these FBVPs with integral boundary conditions and the difficulty associated with finding their exact solutions, some modern numerical algorithms have been developed to find approximate and analytical solutions.

In 2005, Dafterdar-Gejji and Jafari [47] employed the Adomian decomposition method (ADM) to find solutions to a generalized initial system of multi-order fractional differential equations. One year later, they [48] presented an iterative algorithm jointly for solving a general functional equation approximately and called it the Dafterdar-Gejji and Jafari method (DGJIM). Among other numerical algorithms, these two methods, i.e., DGJIM and ADM, are known as two numerical tools with high accuracy and rapid convergence to an exact solution. For more details, one can point out to some works in this regard [49,50,51]. We apply these two strong numerical tools to approximate possible solutions of our suggested FBVP.

In precise terms and with the help of the above ideas, in this paper, we propose a double-order integral FBVP of the multi-term multi-order differential equation in the framework of the Riemann–Liouville (RL) asymmetric derivation operators displayed as

$$\left\{\begin{array}{c}{\displaystyle {\mathfrak{D}}_{0^{+}}^{\varrho }u\left(z\right)=\widehat{\hslash }\left(z,u\left(z\right),{\mathfrak{D}}_{0^{+}}^{{\sigma }_1}u\left(z\right),{\mathfrak{D}}_{0^{+}}^{{\sigma }_2}u\left(z\right),\dots ,{\mathfrak{D}}_{0^{+}}^{{\sigma }_{n-1}}u\left(z\right),{\mathfrak{D}}_{0^{+}}^{{\sigma }_n}u\left(z\right)\right),}\hfill \\ u\left(0\right)=0,u\left(1\right)=p{\mathfrak{I}}_{0^{+}}^{\mu }{k}_1(\xi ,u\left(\xi \right))+q{\mathfrak{I}}_{0^{+}}^{\nu }{k}_2(\eta ,u\left(\eta \right)),\hfill \end{array}\right.$$

(1)

where $0\le z\le 1$, $1<\varrho <2,$$0<{\sigma }_1<{\sigma }_2<\dots <{\sigma }_n<1$, $\varrho >{\sigma }_n+1$, $\widehat{\hslash }:[0,1]\times {\mathbb{R}}^{n+1}\to \mathbb{R}$, $k_j:[0,1]\times \mathbb{R}\to \mathbb{R}$, $(j=1,2)$ are continuous functions; ${\mathfrak{D}}_{0^{+}}^{\varrho },{\mathfrak{D}}_{0^{+}}^{{\sigma }_1},\dots ,{\mathfrak{D}}_{0^{+}}^{{\sigma }_n}$ are RL-derivatives of order $\varrho ,{\sigma }_1,\dots ,{\sigma }_n$, respectively; and ${\mathfrak{I}}_{0^{+}}^{\gamma }$ denotes the RL-integral of order $\gamma \in \{\mu ,\nu \}$ with $\mu ,\nu ,p,q>0$ and $0<\xi ,\eta <1$. Here, we first obtain the corresponding integral equation of the given multi-term multi-order RLFBVP (1) based on a theoretical argument and then establish the existence and uniqueness results with the aid of the fixed point tool. After that, we propose two numerical algorithms entitled DGJIM along with ADM to find approximate solutions.

Indeed, we must emphasize that the novelty and motivation of our work is that, although other papers use the ADM and DGJIM methods for solving IVPs, we here intend to compute approximate solutions for a complicated multi-order multi-term RLFBVP with boundary conditions including double-order RL-fractional integrals. In addition, note that, in the second boundary condition, the value of the unknown function at the end point $z=1$ is proportional to a linear combination of RL-integrals with different orders $\mu ,\nu >0$ at the intermediate points $z=\xi ,\eta \in (0,1)$, respectively. Along with this, we consider the right-hand side nonlinear term $\widehat{\hslash }$ as a multi-variable function including multi-order RL-derivatives finitely.

The rest of this paper is organized as follows. Section 2 recalls fundamental notions on fractional calculus. Section 3 is devoted to establishing some criteria for confirming the existence of solutions. Section 4 introduces the two numerical methods named ADM and DGJIM. In Section 5, the proposed approximation techniques are described using different examples. Some concluding remarks are provided in Section 6.

2. Basic Concepts

First, for the convenience of the readers, we need some fundamental properties and lemmas on fractional calculus which are used further in this paper.

Definition 1.

[3] Let $\varrho >0$ and $\varphi :[0,\infty )\to \mathbb{R}$ be a continuous function. The following integral

$$\left({\mathfrak{I}}_{0^{+}}^{\varrho }\varphi \right)\left(z\right)=\frac{1}{\Gamma \left(\varrho \right)}{\int }_0^z{(z-s)}^{\varrho -1}\varphi \left(s\right)ds,$$

is called the Riemann–Liouville integral of order ϱ such that the integral on the right-hand side exists.

Definition 2.

[3] Let $n-1<\varrho <n$. Then, the ϱth Riemann–Liouville derivative of a continuous function $\varphi :[0,\infty )\to \mathbb{R}$ is defined as

$$\begin{array}{ccc}\hfill {\mathfrak{D}}_{0^{+}}^{\varrho }\varphi \left(z\right)& =& \frac{1}{\Gamma (n-\varrho )}{\left(\frac{d}{dz}\right)}^n{\int }_0^z{(z-s)}^{n-\varrho -1}\varphi \left(s\right)ds\hfill \\ & =& {\left(\frac{d}{dz}\right)}^n{\mathfrak{I}}_{0^{+}}^{n-\varrho }\varphi \left(z\right),\hfill \end{array}$$

provided that the integral on the right-hand side exists and $n=\left[\varrho \right]+1$, where $\left[\varrho \right]$ denotes the greatest integer less than ϱ.

The following properties of the fractional operators are necessary for our paper.

Lemma 1.

[2] Let $u\in L^1(0,1)$ and $\sigma >\varrho >0$. Then,

• ${\displaystyle {\mathfrak{I}}_{0^{+}}^{\sigma }{\mathfrak{I}}_{0^{+}}^{\varrho }u\left(z\right)={\mathfrak{I}}_{0^{+}}^{\sigma +\varrho }u\left(z\right)}$,
• ${\displaystyle {\mathfrak{D}}_{0^{+}}^{\varrho }{\mathfrak{I}}_{0^{+}}^{\sigma }u\left(z\right)={\mathfrak{I}}_{0^{+}}^{\sigma -\varrho }u\left(z\right)}$,
• ${\displaystyle {\mathfrak{D}}_{0^{+}}^{\sigma }{\mathfrak{I}}_{0^{+}}^{\sigma }u\left(z\right)=u\left(z\right)}$.

Lemma 2.

[2] If $\varrho >0$ and $\nu >0$, then

• ${\displaystyle {\mathfrak{D}}_{0^{+}}^{\varrho }{z}^{\nu -1}=\left\{\begin{array}{c}{\displaystyle \frac{\Gamma \left(\nu \right)}{\Gamma (\nu -\varrho )}}{z}^{\nu -\varrho -1},\hfill \\ {\displaystyle {\mathfrak{D}}_{0^{+}}^{\varrho }{z}^{\nu -1}=0,\mathrm{if}\nu -\varrho \in \left\{0\right\}\cup {\mathbb{Z}}^{-},}\hfill \end{array}\right.}$
• ${\displaystyle {\mathfrak{I}}_{0^{+}}^{\varrho }{z}^{\nu }=\frac{\Gamma (\nu +1)}{\Gamma (\nu +\varrho +1)}{z}^{\nu +\varrho }}$.

Lemma 3.

[2] Let $n-1<\varrho <n$ and $u\in C(0,1)$ and ${\mathfrak{D}}_{0^{+}}^{\varrho }u\in L^1(0,1)$. Then,

$${\mathfrak{I}}_{0^{+}}^{\varrho }{\mathfrak{D}}_{0^{+}}^{\varrho }u\left(z\right)=u\left(z\right)-\sum _{j=1}^n\frac{{\mathfrak{I}}_{0^{+}}^{n-\varrho }u\left(0\right)}{\Gamma (\varrho -j+1)}{z}^{\varrho -j},$$

where $n=\left[\varrho \right]+1$ and $\left[\varrho \right]$ denotes the greatest integer less than ϱ.

3. Results of the Existence Criterion

In this section, we first derive an integral equation corresponding to the given multi-term multi-order RLFBVP (1) and then establish required conditions to confirm the existence of solutions for (1).

Definition 3.

The function $u\left(z\right)$ is called a solution for the suggested multi-term multi-order RLFBVP (1) if u satisfies (1) and ${\mathfrak{D}}_{0^{+}}^{\varrho }u\left(z\right)\in C[0,1]$ and $u\left(z\right)\in C[0,1]$.

Theorem 1.

Let $1<\varrho <2,$ $0<{\sigma }_1<{\sigma }_2<\dots <{\sigma }_n<1$, $\varrho >{\sigma }_n+1$, $\mu ,\nu ,p,q>0$, and $0<\xi ,\eta <1$. Then, the function $u\left(z\right)$ is a solution of the RLFBVP (1) if and only if $m\left(z\right)={\mathfrak{D}}_{0^{+}}^{{\sigma }_n}u\left(z\right)$ satisfies the integral equation

$$\begin{array}{cc}\hfill m\left(z\right)& ={\mathfrak{I}}_{0^{+}}^{\varrho -{\sigma }_n}\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right)\hfill \\ \hfill & +\frac{\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)}[\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & +\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & -{\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}]z^{\varrho -{\sigma }_n-1}.\hfill \end{array}$$

(2)

Proof. 

In the first step, let $u\left(z\right)\in C[0,1]$ be a solution of the multi-term multi-order RLFBVP (1) which it gives $m\left(z\right)={\mathfrak{D}}_{0^{+}}^{{\sigma }_n}u\left(z\right)\in C[0,1]$. Applying the RL-operator ${\mathfrak{I}}_{0^{+}}^{{\sigma }_n}$ on both sides of equation $m\left(z\right)={\mathfrak{D}}_{0^{+}}^{{\sigma }_n}u\left(z\right)$, we get

$${\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{\mathfrak{D}}_{0^{+}}^{{\sigma }_n}u\left(z\right)=u\left(z\right)-\frac{\left({\mathfrak{I}}_{0^{+}}^{1-{\sigma }_n}u\right)\left(0\right)}{\Gamma \left({\sigma }_n\right)}{z}^{{\sigma }_n-1}.$$

(3)

Since $\left({\mathfrak{I}}_{0^{+}}^{1-{\sigma }_n}u\right)\left(0\right)=0$, then we have

$$u\left(z\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right).$$

(4)

In view of the second property in Lemma 1 and by (4), it follows that

$$\begin{array}{cc}\hfill {\mathfrak{D}}_{0^{+}}^{{\sigma }_{n-1}}u\left(z\right)& ={\mathfrak{D}}_{0^{+}}^{{\sigma }_{n-1}}{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),\hfill \\ \hfill ⋮& =⋮\hfill \\ \hfill {\mathfrak{D}}_{0^{+}}^{{\sigma }_1}u\left(z\right)& ={\mathfrak{D}}_{0^{+}}^{{\sigma }_1}{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right).\hfill \end{array}$$

Since $1<\varrho <2$, by definition of the Riemann–Liouville fractional derivative, ${\mathfrak{D}}_{0^{+}}^{\varrho }u\left(z\right)={\mathfrak{D}}_{0^{+}}^2{\mathfrak{I}}_{0^{+}}^{2-\varrho }u\left(z\right)$. Now, by (4), we get ${\mathfrak{D}}_{0^{+}}^{\varrho }u\left(z\right)={\mathfrak{D}}_{0^{+}}^2{\mathfrak{I}}_{0^{+}}^{2-\varrho }{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right)$. Now, by Lemma 1, if we use the semi-group property for Riemann–Liouville fractional integrals, we have

$${\mathfrak{I}}_{0^{+}}^{2-\varrho }{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right)={\mathfrak{I}}_{0^{+}}^{2-\varrho +{\sigma }_n}m\left(z\right).$$

Again, by definition of the Riemann–Liouville fractional derivative, we have

$${\mathfrak{D}}_{0^{+}}^2{\mathfrak{I}}_{0^{+}}^{2-\varrho +{\sigma }_n}m\left(z\right)={\mathfrak{D}}_{0^{+}}^{-(-\varrho +{\sigma }_n)}m\left(z\right)={\mathfrak{D}}_{0^{+}}^{\varrho -{\sigma }_n}m\left(z\right),$$

and so

$${\mathfrak{D}}_{0^{+}}^{\varrho }u\left(z\right)={\mathfrak{D}}_{0^{+}}^{\varrho -{\sigma }_n}m\left(z\right).$$

Consequently, the multi-term multi-order equation illustrated by (1), becomes

$${\mathfrak{D}}_{0^{+}}^{\varrho -{\sigma }_n}m\left(z\right)=\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right),0\le z\le 1.$$

(5)

Setting $\lambda =\varrho -{\sigma }_n>1$, ${\lambda }_j={\sigma }_n-{\sigma }_j$, ${\sigma }_0=0$ $(j=0,1,\dots n)$, then (5) can be rewritten as

$${\mathfrak{D}}_{0^{+}}^{\lambda }m\left(z\right)=\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\lambda }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}m\left(z\right),m\left(z\right)\right),0\le z\le 1.$$

(6)

Hence, by (4), it follows that $u\left(0\right)=0$, and one can determine the value of the initial condition $m\left(0\right)$. Therefore, since $m\left(z\right)\in C[0,1]$,

$${\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right)=\frac{1}{\Gamma \left({\sigma }_n\right)}{\int }_0^z{(z-s)}^{{\sigma }_n-1}m\left(s\right)ds,$$

and so we can arbitrarily provide the initial value of $m\left(z\right)$ such that $u\left(0\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right){|}_{z=0}=0$. We assume that

$$m\left(0\right)=0.$$

(7)

Now, taking the Riemann–Liouville fractional integral ${\mathfrak{I}}_{0^{+}}^{\lambda }$ on both sides of (6), we find that

$${\mathfrak{I}}_{0^{+}}^{\lambda }{\mathfrak{D}}_{0^{+}}^{\lambda }m\left(z\right)={\mathfrak{I}}_{0^{+}}^{\lambda }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\lambda }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}m\left(z\right),m\left(z\right)\right),0\le z\le 1.$$

(8)

By the hypothesis of the theorem, we have $\lambda =\varrho -{\sigma }_n>1$. Then, from Lemma 3, the left-hand side of (8) becomes

$${\mathfrak{I}}_{0^{+}}^{\lambda }{\mathfrak{D}}_{0^{+}}^{\lambda }m\left(z\right)=m\left(z\right)+c_1{z}^{\lambda -1}+c_2{z}^{\lambda -2},$$

hence Equation (8) is rewritten in the following form

$$m\left(z\right)={\mathfrak{I}}_{0^{+}}^{\lambda }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\lambda }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}m\left(z\right),m\left(z\right)\right)-c_1{z}^{\lambda -1}-c_2{z}^{\lambda -2}.$$

(9)

By (7), since $m\left(0\right)=0$ and $2>\lambda >1$, we get $c_2=0$. Therefore, Equation (9) becomes

$$m\left(z\right)={\mathfrak{I}}_{0^{+}}^{\lambda }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\lambda }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}m\left(z\right),m\left(z\right)\right)-c_1{z}^{\lambda -1}.$$

(10)

By using the second boundary condition given in (1) and by (4), we have

$$\begin{array}{c}\hfill u\left(1\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right){|}_{z=1}=p{\mathfrak{I}}_{0^{+}}^{\mu }{k}_1(\xi ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(\xi \right))+q{\mathfrak{I}}_{0^{+}}^{\nu }{k}_2(\eta ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(\eta \right)).\end{array}$$

(11)

With the help of Lemma 1 and from (10) and (11), we figure out that

$$\begin{array}{cc}\hfill u\left(1\right)& ={\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right){|}_{z=1}\hfill \\ \hfill & ={\mathfrak{I}}_{0^{+}}^{{\sigma }_n+\lambda }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\lambda }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}-c_1{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{z}^{\lambda -1}{|}_{z=1}\hfill \\ \hfill & ={\mathfrak{I}}_{0^{+}}^{{\sigma }_n+\lambda }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\lambda }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}-c_1\frac{\Gamma \left(\lambda \right)}{\Gamma (\lambda +{\sigma }_n)}{z}^{\lambda +{\sigma }_n-1}{|}_{z=1}\hfill \\ \hfill & =\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right))ds+\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right))ds.\hfill \end{array}$$

However, we have $\lambda +{\sigma }_n-1=\varrho -{\sigma }_n+{\sigma }_n-1=\varrho -1>0$. Then, one can write

$$\begin{array}{cc}\hfill & \frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right))ds+\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right))ds\hfill \\ \hfill & ={\mathfrak{I}}_{0^{+}}^{{\sigma }_n+\lambda }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\lambda }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}-c_1\frac{\Gamma \left(\lambda \right)}{\Gamma (\lambda +{\sigma }_n)}\hfill \\ \hfill & ={\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}-c_1\frac{\Gamma (\varrho -{\sigma }_n)}{\Gamma \left(\varrho \right)}.\hfill \end{array}$$

Thus, we get

$$\begin{array}{cc}\hfill c_1& =\frac{\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)}[{\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}\hfill \\ \hfill & -\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right))ds\hfill \\ \hfill & -\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right))ds].\hfill \end{array}$$

By substituting the value of $c_1$ into Equation (10), we obtain the following equation

$$\begin{array}{cc}\hfill m\left(z\right)& ={\mathfrak{I}}_{0^{+}}^{\varrho -{\sigma }_n}\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right)\hfill \\ \hfill & +\frac{\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)}[\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & +\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & -{\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}]z^{\varrho -{\sigma }_n-1},\hfill \end{array}$$

which implies that $m\left(z\right)={\mathfrak{D}}_{0^{+}}^{{\sigma }_n}u\left(z\right)\in C[0,1]$ is a solution of (2).

Conversely, suppose that $m\left(z\right)={\mathfrak{D}}_{0^{+}}^{{\sigma }_n}u\left(z\right)\in C[0,1]$ is a solution of (2). By applying the Riemann–Liouville fractional integral ${\mathfrak{I}}_{0^{+}}^{{\sigma }_n}$ on both sides of $m\left(z\right)={\mathfrak{D}}_{0^{+}}^{{\sigma }_n}u\left(z\right)$, we have

$${\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{\mathfrak{D}}_{0^{+}}^{{\sigma }_n}u\left(z\right)=u\left(z\right)-\frac{\left({\mathfrak{I}}_{0^{+}}^{1-{\sigma }_n}u\right)\left(0\right)}{\Gamma \left({\sigma }_n\right)}{z}^{{\sigma }_n-1}.$$

Due to $\left({\mathfrak{I}}_{0^{+}}^{1-{\sigma }_n}u\right)\left(0\right)=0$, we obtain $u\left(z\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right).$ In the next steps, we obtain other fractional derivatives recursively and the second property in Lemma 1 as follows

$$\begin{array}{cc}\hfill u\left(z\right)& ={\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),\hfill \\ \hfill {\mathfrak{D}}_{0^{+}}^{{\sigma }_{n-1}}u\left(z\right)& ={\mathfrak{D}}_{0^{+}}^{{\sigma }_{n-1}}{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),\hfill \\ \hfill ⋮& =⋮\hfill \\ \hfill {\mathfrak{D}}_{0^{+}}^{{\sigma }_1}u\left(z\right)& ={\mathfrak{D}}_{0^{+}}^{{\sigma }_1}{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right).\hfill \end{array}$$

(12)

By taking the Riemann–Liouville operator ${\mathfrak{I}}_{0^{+}}^{{\sigma }_n}$ on both sides of (2), it becomes

$$\begin{array}{cc}\hfill {\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right)& ={\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right)\hfill \\ \hfill & +\frac{\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)}[\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & +\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & -{\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}]{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{z}^{\varrho -{\sigma }_n-1},\hfill \end{array}$$

and so

$$\begin{array}{cc}\hfill u\left(z\right)& ={\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right)\hfill \\ \hfill & +\frac{\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)}[\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & +\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & -{\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}]{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{z}^{\varrho -{\sigma }_n-1}.\hfill \end{array}$$

(13)

In the sequel, by applying the Riemann–Liouville operator ${\mathfrak{D}}_{0^{+}}^{\varrho }$ on both sides of (13), it follows

$$\begin{array}{cc}\hfill {\mathfrak{D}}_{0^{+}}^{\varrho }u\left(z\right)& ={\mathfrak{D}}_{0^{+}}^{\varrho }{\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right)\hfill \\ \hfill & +\frac{\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)}[\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & +\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & -{\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}]{\mathfrak{D}}_{0^{+}}^{\varrho }{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{z}^{\varrho -{\sigma }_n-1}.\hfill \end{array}$$

Since, by Lemma 2, ${\displaystyle {\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{z}^{\varrho -{\sigma }_n-1}=\frac{\Gamma (\varrho -{\sigma }_n)}{\Gamma \left(\varrho \right)}{z}^{\varrho -1}}$ and ${\mathfrak{D}}_{0^{+}}^{\varrho }{z}^{\varrho -1}=0$, we get

$$\begin{array}{cc}\hfill {\mathfrak{D}}_{0^{+}}^{\varrho }u\left(z\right)& =\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right)\hfill \\ \hfill & +[\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & +\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & -{\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}]{\mathfrak{D}}_{0^{+}}^{\varrho }{z}^{\varrho -1}\hfill \\ \hfill & =\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right).\hfill \end{array}$$

(14)

According to (12), the fractional differential Equation (14) reduces to

$${\mathfrak{D}}_{0^{+}}^{\varrho }u\left(z\right)=\widehat{\hslash }\left(z,u\left(z\right),{\mathfrak{D}}_{0^{+}}^{{\sigma }_1}u\left(z\right),{\mathfrak{D}}_{0^{+}}^{{\sigma }_2}u\left(z\right),\dots ,{\mathfrak{D}}_{0^{+}}^{{\sigma }_{n-1}}u\left(z\right),{\mathfrak{D}}_{0^{+}}^{{\sigma }_n}u\left(z\right)\right).$$

Finally, we check both boundary conditions of problem (1). In view of Equation (2) and by definition of the Riemann–Liouville integral of the function

$$\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right)$$

of order $\varrho -{\sigma }_n$ at point $z=0$, it is immediately deduced that

$$\begin{array}{cc}\hfill m\left(0\right)& ={\mathfrak{I}}_{0^{+}}^{\varrho -{\sigma }_n}\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=0}\hfill \\ \hfill & +\frac{\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)}[\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & +\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & -{\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}]z^{\varrho -{\sigma }_n-1}{|}_{z=0}\hfill \\ \hfill & =0+0=0.\hfill \end{array}$$

(15)

Thus, $m\left(0\right)=0$. Hence, we have $u\left(z\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right)$, and so $u\left(0\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right){|}_{z=0}=0$. Thus, $u\left(0\right)=0$. This means that the first boundary condition holds. Now, to check the second boundary condition, by substituting $z=1$ into (13), we obtain

$$\begin{array}{cc}\hfill u\left(1\right)& ={\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}\hfill \\ \hfill & +\frac{\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)}[\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & +\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & -{\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}]{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{z}^{\varrho -{\sigma }_n-1}{|}_{z=1}\hfill \\ \hfill & ={\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}\hfill \\ \hfill & +\frac{\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)}[\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & +\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & -{\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}]\frac{\Gamma (\varrho -{\sigma }_n)}{\Gamma \left(\varrho \right)}{z}^{\varrho -1}{|}_{z=1}\hfill \\ \hfill & =p{\mathfrak{I}}_{0^{+}}^{\mu }{k}_1(\xi ,u\left(\xi \right))+q{\mathfrak{I}}_{0^{+}}^{\nu }{k}_2(\eta ,u\left(\eta \right)).\hfill \end{array}$$

Therefore, we figure out that $u\left(z\right)$ satisfies the multi-term multi-order RLFBVP (1) and so u will be a solution of the mentioned RLFBVP, and the proof is completed. □

Here, we introduce the Banach space $E=C[0,1]$ with the norm ${\displaystyle \parallel m\parallel =\underset{z\in [0,1]}{max}\left|m\left(z\right)\right|}$, and, along with this, by Theorem 1, we define an operator $\Psi :E⟶E$ by

$$\begin{array}{cc}\hfill (\Psi m)\left(z\right)& ={\mathfrak{I}}_{0^{+}}^{\varrho -{\sigma }_n}\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right)\hfill \\ \hfill & +\frac{\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)}[\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & +\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(s\right)\right)ds\hfill \\ \hfill & -{\mathfrak{I}}_{0^{+}}^{\varrho }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right),{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_1}m\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n-{\sigma }_{n-1}}m\left(z\right),m\left(z\right)\right){|}_{z=1}]z^{\varrho -{\sigma }_n-1}.\hfill \end{array}$$

(16)

We clearly have the following equation

$$\Psi m=m,m\in E,$$

(17)

which is equivalent to Equation (2). If $\Psi$ has a fixed point, then it will be the solution of the multi-term multi-order RLFBVP (1). On the other side, notice that the continuity of all three functions $\widehat{\hslash },k_1$, and $k_2$ confirms that of the operator $\Psi$. In this place, we want to express the existence theorem in relation to solutions of the multi-term multi-order RLFBVP (1).

Theorem 2.

Assume that these assumptions are valid:

• ($\mathcal{AS}1$) There exist real constants $M_j(j=0,1,\dots ,n)$ such that

  $$|\widehat{\hslash }(z,u_0,u_1,\dots ,u_n)-\widehat{\hslash }(z,U_0,U_1,\dots ,U_n)|\le \sum _{j=0}^n{M}_j|u_j-U_j|,$$

  for all $z\in [0,1]$ and $(u_0,u_1,\dots ,u_n),(U_0,U_1,\dots ,U_n)\in {\mathbb{R}}^{n+1}$.
• ($\mathcal{AS}2$) There exist two real constants ${\theta }_1,{\theta }_2>0$ such that

  $$\begin{array}{cc}\hfill |k_1(z,m)-k_1(z,u)|& \le {\theta }_1|m-u|,m,u\in \mathbb{R},\hfill \\ \hfill |k_2(z,m)-k_2(z,u)|& \le {\theta }_2|m-u|,m,u\in \mathbb{R}.\hfill \end{array}$$
• ($\mathcal{AS}3$) Let

  $$\begin{array}{cc}\hfill 0<\Phi & =\frac{\Gamma \left(\varrho \right)p{\theta }_1{\xi }^{\mu }}{\Gamma (\varrho -{\sigma }_n)\Gamma (\mu +1)\Gamma ({\sigma }_n+1)}+\frac{\Gamma \left(\varrho \right)q{\theta }_2{\eta }^{\nu }}{\Gamma (\varrho -{\sigma }_n)\Gamma (\nu +1)\Gamma ({\sigma }_n+1)}\hfill \\ \hfill & +\sum _{j=0}^n\left[\frac{M_j}{\Gamma (\varrho -{\sigma }_j+1)}+\frac{M_j\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)\Gamma (\varrho +{\sigma }_n-{\sigma }_j+1)}\right]<1.\hfill \end{array}$$

Then, the multi-term multi-order RLFBVP (1) has a unique solution.

Proof. 

In view of Theorem 1, it is explicit that the existence of solutions to the multi-term multi-order RLFBVP (1) is derived from the existence of solutions to Equation (16) or (17). Thus, it suffices to prove that (16) has a unique fixed point. Now, let $\lambda =\varrho -{\sigma }_n$, ${\sigma }_0=0$, and ${\lambda }_j={\sigma }_n-{\sigma }_j$ for $j=0,1,\dots ,n$. Then, from $\left(\mathcal{AS}1\right)$, it follows that for any $m_1,m_2\in E,$ we have

$$\begin{array}{cc}\hfill & |\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}{m}_1\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}{m}_1\left(z\right),m_1\left(z\right)\right)-\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}{m}_2\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}{m}_2\left(z\right),m_2\left(z\right)\right)|\hfill \\ \hfill & \le \sum _{j=0}^n{M}_j|{\mathfrak{I}}_{0^{+}}^{{\lambda }_j}{m}_1\left(z\right)-{\mathfrak{I}}_{0^{+}}^{{\lambda }_j}{m}_2\left(z\right)|.\hfill \end{array}$$

(18)

Taking the Riemann–Liouville operator ${\mathfrak{I}}_{0^{+}}^{\lambda }$ on both sides of inequality (18), we find that

$$\begin{array}{cc}\hfill {\mathfrak{I}}_{0^{+}}^{\lambda }\left|\widehat{\hslash }\right(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}{m}_1\left(z\right),\dots ,& {\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}{m}_1\left(z\right),m_1\left(z\right))-\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}{m}_2\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}{m}_2\left(z\right),m_2\left(z\right)\right)|\hfill \\ \hfill & \le {\mathfrak{I}}_{0^{+}}^{\lambda }\sum _{j=0}^n{M}_j|{\mathfrak{I}}_{0^{+}}^{{\lambda }_j}{m}_1\left(z\right)-{\mathfrak{I}}_{0^{+}}^{{\lambda }_j}{m}_2\left(z\right)|\le \sum _{j=0}^n{M}_j{\mathfrak{I}}_{0^{+}}^{\lambda +{\lambda }_j}|m_1\left(z\right)-m_2\left(z\right)|\hfill \\ \hfill & \le \parallel m_1-m_2\parallel \sum _{j=0}^n\frac{M_j}{\Gamma (\lambda +{\lambda }_j+1)}=\parallel m_1-m_2\parallel \sum _{j=0}^n\frac{M_j}{\Gamma (\varrho -{\sigma }_n+{\sigma }_n-{\sigma }_j+1)}\hfill \\ \hfill & =\parallel m_1-m_2\parallel \sum _{j=0}^n\frac{M_j}{\Gamma (\varrho -{\sigma }_j+1)}.\hfill \end{array}$$

(19)

On the other side, by using $\left(\mathcal{AS}2\right)$, we get

$$\begin{array}{cc}\hfill |\frac{\Gamma \left(\varrho \right)}{\Gamma \left(\lambda \right)}& [\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{m}_1\left(s\right)\right)ds+\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{m}_1\left(s\right)\right)ds\hfill \\ \hfill & -{\mathfrak{I}}_{0^{+}}^{{\sigma }_n+\lambda }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}{m}_1\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}{m}_1\left(z\right),m_1\left(z\right)\right){|}_{z=1}]z^{\lambda -1}\hfill \\ \hfill & -\frac{\Gamma \left(\varrho \right)}{\Gamma \left(\lambda \right)}[\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}{k}_1\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{m}_2\left(s\right)\right)ds+\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}{k}_2\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{m}_2\left(s\right)\right)ds\hfill \\ \hfill & -{\mathfrak{I}}_{0^{+}}^{{\sigma }_n+\lambda }\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}{m}_2\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}{m}_2\left(z\right),m_2\left(z\right)\right){|}_{z=1}\left]z^{\lambda -1}\right|\hfill \\ \hfill & \le \frac{\Gamma \left(\varrho \right)}{\Gamma \left(\lambda \right)}\left[\frac{p}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}\right|k_1\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{m}_1\left(s\right)\right)-k_1(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{m}_2\left(s\right)|ds\hfill \\ \hfill & +\frac{q}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}|k_2\left(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{m}_1\left(s\right)\right)-k_2(s,{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{m}_2\left(s\right)|ds\hfill \\ \hfill & +{\mathfrak{I}}_{0^{+}}^{{\sigma }_n+\lambda }|\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}{m}_1\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}{m}_1\left(z\right),m_1\left(z\right)\right)-\widehat{\hslash }\left(z,{\mathfrak{I}}_{0^{+}}^{{\lambda }_0}{m}_2\left(z\right),\dots ,{\mathfrak{I}}_{0^{+}}^{{\lambda }_{n-1}}{m}_2\left(z\right),m_2\left(z\right)\right)|{|}_{z=1}]\hfill \\ \hfill & \le \frac{\Gamma \left(\varrho \right)}{\Gamma \left(\lambda \right)}\left[\frac{p{\theta }_1}{\Gamma \left(\mu \right)}{\int }_0^{\xi }{(\xi -s)}^{\mu -1}\right|{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{m}_1\left(s\right)-{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{m}_2\left(s\right)|ds\hfill \\ \hfill & +\frac{q{\theta }_2}{\Gamma \left(\nu \right)}{\int }_0^{\eta }{(\eta -s)}^{\nu -1}|{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{m}_1\left(s\right)-{\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{m}_2\left(s\right)|ds+{\mathfrak{I}}_{0^{+}}^{{\sigma }_n+\lambda }\sum _{j=0}^n{M}_j{\mathfrak{I}}_{0^{+}}^{{\lambda }_j}|m_1\left(z\right)-m_2\left(z\right)|{|}_{z=1}]\hfill \\ \hfill & \le \frac{\Gamma \left(\varrho \right)}{\Gamma \left(\lambda \right)}\left[\frac{p{\theta }_1{\xi }^{\mu }}{\Gamma (\mu +1)\Gamma ({\sigma }_n+1)}+\frac{q{\theta }_2{\eta }^{\nu }}{\Gamma (\nu +1)\Gamma ({\sigma }_n+1)}+\sum _{j=0}^n\frac{M_j}{\Gamma ({\sigma }_n+\lambda +{\lambda }_j+1)}\right]\parallel m_1-m_2\parallel \hfill \\ \hfill & =\frac{\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)}\left[\frac{p{\theta }_1{\xi }^{\mu }}{\Gamma (\mu +1)\Gamma ({\sigma }_n+1)}+\frac{q{\theta }_2{\eta }^{\nu }}{\Gamma (\nu +1)\Gamma ({\sigma }_n+1)}+\sum _{j=0}^n\frac{M_j}{\Gamma ({\sigma }_n+\varrho -{\sigma }_j+1)}\right]\parallel m_1-m_2\parallel .\hfill \end{array}$$

(20)

Consequently, by adding both sides of (19) and (20) and according to the definition of $\Psi$ in (16), we have

$$\begin{array}{cc}\hfill |\Psi m_1\left(z\right)-\Psi m_2\left(z\right)|& \le [\frac{\Gamma \left(\varrho \right)p{\theta }_1{\xi }^{\mu }}{\Gamma (\varrho -{\sigma }_n)\Gamma (\mu +1)\Gamma ({\sigma }_n+1)}+\frac{\Gamma \left(\varrho \right)q{\theta }_2{\eta }^{\nu }}{\Gamma (\varrho -{\sigma }_n)\Gamma (\nu +1)\Gamma ({\sigma }_n+1)}\hfill \\ \hfill & +\sum _{j=0}^n\left(\frac{M_j}{\Gamma (\varrho -{\sigma }_j+1)}+\frac{M_j\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)\Gamma (\varrho +{\sigma }_n-{\sigma }_j+1)}\right)]\parallel m_1-m_2\parallel .\hfill \end{array}$$

By using $\left(\mathcal{AS}3\right)$, we find

$$\parallel \Psi m_1-\Psi m_2\parallel \le \Phi \parallel m_1-m_2\parallel ,$$

where $\Phi \in (0,1)$. Hence, by the Banach fixed point theorem [52], it follows that $\Psi$ has a unique fixed point which points out that the suggested multi-term multi-order RLFBVP (1) has a unique solution. □

4. Approximation of Solutions via DGJIM and ADM Methods

This section is devoted to implementing the numerical methods named DGJIM and ADM. Indeed, we here state how we can employ these methods to our suggested multi-term multi-order RLFBVP. In both algorithms, appropriate recursion relations are formulated to approximate the solutions of (1) along with their convergence. Our techniques are inspired by [47,48].

4.1. DGJIM Numerical Method

We prove above that the solutions of Equations (1) and (2) are equivalent. Thus, we now suppose that the right-hand side of (17) is written under the following decomposition (not uniquely)

$$(\Psi m)\left(z\right)=\tilde{\mathbf{L}}\left(m\left(z\right)\right)+\tilde{\mathbf{N}}\left(m\left(z\right)\right)+\zeta \left(z\right),$$

where the operator $\tilde{\mathbf{L}}$ is linear, the operator $\tilde{\mathbf{N}}$ stands for the nonlinear terms, and $\zeta$ is a known function. Then, one can rewrite (2) in the decomposed form

$$\begin{array}{c}\hfill m\left(z\right)=\tilde{\mathbf{L}}\left(m\left(z\right)\right)+\tilde{\mathbf{N}}\left(m\left(z\right)\right)+\zeta \left(z\right).\end{array}$$

(21)

Suppose that the solution of (21) is written as a series as follows

$$\begin{array}{c}\hfill m\left(z\right)=\sum _{n=0}^{+\infty }{m}_n\left(z\right).\end{array}$$

(22)

By combining (22) and (21), we get

$$\begin{array}{c}\hfill \sum _{n=0}^{+\infty }{m}_n\left(z\right)=\tilde{\mathbf{L}}\left(\sum _{n=0}^{+\infty }{m}_n\left(z\right)\right)+\tilde{\mathbf{N}}\left(\sum _{n=0}^{+\infty }{m}_n\left(z\right)\right)+\zeta \left(z\right).\end{array}$$

(23)

Since $\tilde{\mathbf{L}}$ is linear, by a simple manipulation, we obtain the following algorithm known as the DGJIM numerical method:

$$\left\{\begin{array}{c}{m}_0\left(z\right)=\zeta \left(z\right),\hfill \\ m_1\left(z\right)=\tilde{\mathbf{L}}\left(m_0\left(z\right)\right)+\tilde{\mathbf{N}}\left(m_0\left(z\right)\right),\hfill \\ m_2\left(z\right)=\tilde{\mathbf{L}}\left(m_1\left(z\right)\right)+\tilde{\mathbf{N}}\left(m_0\left(z\right)+m_1\left(z\right)\right)-\tilde{\mathbf{N}}\left(m_0\left(z\right)\right),\hfill \\ m_3\left(z\right)=\tilde{\mathbf{L}}\left(m_2\left(z\right)\right)+\tilde{\mathbf{N}}\left(m_0\left(z\right)+m_1\left(z\right)+m_2\left(z\right)\right)-\tilde{\mathbf{N}}\left(m_0\left(z\right)+m_1\left(z\right)\right),\hfill \\ ⋮=⋮\hfill \\ {\displaystyle m_n\left(z\right)=\tilde{\mathbf{L}}\left(m_{n-1}\left(z\right)\right)+\tilde{\mathbf{N}}\left(\sum _{i=0}^{n-1}{m}_i\left(z\right)\right)-\tilde{\mathbf{N}}\left(\sum _{i=0}^{n-2}{m}_i\left(z\right)\right),}\hfill \\ ⋮=⋮\hfill \end{array}\right.$$

(24)

Therefore, we can obtain the n-term approximate solution of the integral Equation (2) as

$$w_n\left(z\right)=\sum _{i=0}^n{m}_i\left(z\right).$$

(25)

In view of (25), we simply get

$$\begin{array}{c}\hfill m_n\left(z\right)=w_n\left(z\right)-w_{n-1}\left(z\right).\end{array}$$

(26)

Thus, a combination of (24) and (26) gives

$$w_n\left(z\right)=w_{n-1}\left(z\right)+\tilde{\mathbf{L}}\left(w_{n-1}\left(z\right)-w_{n-2}\left(z\right)\right)+\tilde{\mathbf{N}}\left(w_{n-1}\left(z\right)\right)-\tilde{\mathbf{N}}\left(w_{n-2}\left(z\right)\right).$$

(27)

Now, let

$$\begin{array}{cc}\hfill \parallel \tilde{\mathbf{L}}m-\tilde{\mathbf{L}}u\parallel & \le {\mu }_1\parallel m-u\parallel ,0<{\mu }_1<1,\hfill \\ \hfill \parallel \tilde{\mathbf{N}}m-\tilde{\mathbf{N}}u\parallel & \le {\mu }_2\parallel m-u\parallel ,0<{\mu }_2<1,\hfill \end{array}$$

where ${\mu }_1+{\mu }_2<1.$ Therefore, the Banach fixed point principle guarantees the existence of a unique solution $\tilde{w}\left(z\right)$ for (21) and so for the integral Equation (2). According to the relation (27), the following iterative expression is derived

$$\begin{array}{cc}\hfill \parallel w_n-w_{n-1}\parallel & \le {\mu }_1\parallel w_{n-1}-w_{n-2}\parallel +{\mu }_2\parallel w_{n-1}-w_{n-2}\parallel \hfill \\ \hfill & =({\mu }_1+{\mu }_2)\parallel w_{n-1}-w_{n-2}\parallel \hfill \\ \hfill & \le {({\mu }_1+{\mu }_2)}^2\parallel w_{n-2}-w_{n-3}\parallel \hfill \\ \hfill & \le ⋮⋮⋮\hfill \\ \hfill & \le {\left({\mu }_1+{\mu }_2\right)}^{n-1}\parallel w_1-w_0\parallel ,\hfill \end{array}$$

which implies the absolute convergence and the uniform convergence of the sequence $\left\{w_n\right\}$ to the exact solution $\tilde{w}\left(z\right)$.

4.2. ADM Numerical Method

To implement the ADM numerical method, the nonlinear term ${\displaystyle \tilde{\mathbf{N}}\left(\sum _{n=0}^{+\infty }{m}_n\left(z\right)\right)}$ introduced in (23) is decomposed into a series of Adomian polynomials as

$$\tilde{\mathbf{N}}\left(\sum _{n=0}^{+\infty }{m}_n\left(z\right)\right)=\sum _{n=0}^{+\infty }{\mathbf{A}}_n\left(m_0,m_1,\dots ,m_n\right),$$

where ${\mathbf{A}}_n\left(m_0,m_1,\dots ,m_n\right)$ is produced by

$${\mathbf{A}}_n\left(m_0,m_1,\dots ,m_n\right)=\frac{1}{n!}\frac{{\partial }^n}{\partial z^n}{\left[\tilde{\mathbf{N}}\left(\sum _{k=0}^{+\infty }{m}_k{z}^k\right)\right]}_{z=0},(n\in \mathbb{N}\cup \left\{0\right\}).$$

(28)

Consequently, Equation (23) reduces to

$$\sum _{n=0}^{+\infty }{m}_n\left(z\right)=\tilde{\mathbf{L}}\left(\sum _{n=0}^{+\infty }{m}_n\left(z\right)\right)+\sum _{n=0}^{+\infty }{\mathbf{A}}_n\left(m_0\left(z\right),m_1\left(z\right),\dots ,m_n\left(z\right)\right)+\zeta \left(z\right),$$

which gives us the following iterative schemes called the ADM method:

$$\left\{\begin{array}{c}{m}_0\left(z\right)=\zeta \left(z\right),\hfill \\ m_1\left(z\right)=\tilde{\mathbf{L}}\left(m_0\left(z\right)\right)+{\mathbf{A}}_0\left(m_0\left(z\right),m_1\left(z\right),\dots ,m_n\left(z\right)\right),\hfill \\ m_2\left(z\right)=\tilde{\mathbf{L}}\left(m_1\left(z\right)\right)+{\mathbf{A}}_1\left(m_0\left(z\right),m_1\left(z\right),\dots ,m_n\left(z\right)\right),\hfill \\ m_3\left(z\right)=\tilde{\mathbf{L}}\left(m_2\left(z\right)\right)+{\mathbf{A}}_2\left(m_0\left(z\right),m_1\left(z\right),\dots ,m_n\left(z\right)\right),\hfill \\ ⋮⋮⋮⋮⋮⋮\hfill \\ m_n\left(z\right)=\tilde{\mathbf{L}}\left(m_{n-1}\left(z\right)\right)+{\mathbf{A}}_{n-1}\left(m_0\left(z\right),m_1\left(z\right),\dots ,m_n\left(z\right)\right),\hfill \\ ⋮⋮⋮⋮⋮⋮\hfill \end{array}\right.$$

(29)

Finally, by writing M-term approximate solution of the integral Equation (2) as

$$w_M\left(z\right)=\sum _{n=0}^M{m}_n\left(z\right),$$

(30)

we obtain the exact solution of (2) by

$$m\left(z\right)=\underset{M\to +\infty }{lim}{w}_M\left(z\right).$$

(31)

Lastly, we find that the approximate solutions and the exact solution of the multi-term multi-order RLFBVP (1) are extracted as $u_n\left(z\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n}{w}_n\left(z\right)$ and $u\left(z\right)={\mathfrak{I}}_{0^{+}}^{{\sigma }_n}m\left(z\right)$, respectively.

5. Application

Here, we prepare two distinct examples. In the first, the theoretical existence results are examined, and, in the second, the approximate solutions of a given RLFBVP are obtained with the help of the DGJIM and ADM numerical methods introduced above. Note that, in the second example, we compare the approximate solutions obtained by two mentioned numerical methods with the exact ones for different given fractional orders.

Example 1.

Let us consider the following RLFBVP

$$\left\{\begin{array}{c}{\displaystyle {\mathfrak{D}}_{0^{+}}^{1.8}u\left(z\right)=z^2+\frac{1}{8}sin\left(2u\left(z\right)\right)+\frac{1}{4}{\mathfrak{D}}_{0^{+}}^{0.4}u\left(z\right)+\frac{2}{10}arctan\left({\mathfrak{D}}_{0^{+}}^{0.5}u\left(z\right)\right),z\in (0,1),}\hfill \\ u\left(0\right)=0,\hfill \\ {\displaystyle u\left(1\right)=6{\int }_0^{\frac{1}{2}}\frac{{(1-2s)}^3\left(1+u\left(s\right)\right)}{8\Gamma \left(4\right)(4+s^2)}ds+24{\int }_0^{\frac{1}{4}}\frac{{(1-4s)}^4\left(e^{-s}+sin\left(u\left(s\right)\right)\right)}{\Gamma \left(5\right)1024}ds,}\hfill \end{array}\right.$$

where we take data $\varrho =1.8$, $n=2$, ${\sigma }_0=0$, ${\sigma }_1=0.4$, ${\sigma }_2=0.5$, $\xi ={\displaystyle \frac{1}{2}}$, $\eta ={\displaystyle \frac{1}{4}}$, $p=6$, $q=24$, $\mu =4$, and $\nu =5$. Along with these, continuous functions

$$\hslash (z,s\left(z\right),x\left(z\right),y\left(z\right))=z^2+\frac{1}{8}sin\left(2s\left(z\right)\right)+\frac{1}{4}x\left(z\right)+\frac{2}{10}arctan\left(y\left(z\right)\right),$$

and

$$k_1(z,u\left(z\right))=\frac{1+u\left(z\right)}{4+z^2},k_2(z,u\left(z\right))=\frac{e^{-z}+sin\left(u\left(z\right)\right)}{4},$$

are defined on their domain. Clearly, $M_0=M_1=0.25$ and $M_2=0.2$. On the other side, we get

$$|k_1(z,u\left(z\right))-k_1(z,U\left(z\right))|\le |\frac{1+u\left(z\right)}{4+z^2}-\frac{1+U\left(z\right)}{4+z^2}|\le \frac{1}{4+z^2}|u\left(z\right)-U\left(z\right)|,$$

and

$$|k_2(z,u\left(z\right))-k_2(z,U\left(z\right))|\le |\frac{e^{-z}+sin\left(u\left(z\right)\right)}{4}-\frac{e^{-z}+sin\left(U\left(z\right)\right)}{4}|\le \frac{1}{4}|u\left(z\right)-U\left(z\right)|.$$

Thus, ${\theta }_1={\theta }_2=0.25$. In addition,

$$\begin{array}{cc}\hfill \Phi & =\frac{\Gamma \left(\varrho \right)p{\theta }_1{\xi }^{\mu }}{\Gamma (\varrho -{\sigma }_n)\Gamma (\mu +1)\Gamma ({\sigma }_n+1)}+\frac{\Gamma \left(\varrho \right)q{\theta }_2{\eta }^{\nu }}{\Gamma (\varrho -{\sigma }_n)\Gamma (\nu +1)\Gamma ({\sigma }_n+1)}\hfill \\ \hfill & +\sum _{j=0}^n\left[\frac{M_j}{\Gamma (\varrho -{\sigma }_j+1)}+\frac{M_j\Gamma \left(\varrho \right)}{\Gamma (\varrho -{\sigma }_n)\Gamma (\varrho +{\sigma }_n-{\sigma }_j+1)}\right]\approx 0.8951<1\hfill \end{array}$$

In consequence, by Theorem 2, a unique solution exists for the multi-term multi-order RLFBVP considered above.

For the next example, we consider three different cases for the order of the proposed RLFBVP and compare obtained approximate results with exact outcomes, which shows the effectiveness of both DGJIM and ADM numerical methods together.

Example 2.

In the present example, we consider three distinct values for ϱ as $\varrho =1.4$, $\varrho =1.7$ and $\varrho =1.9$.

$•\mathbf{Case}\left(\mathbf{I}\right):\varrho =1.4$: Let us consider the following RLFBVP which has a structure as

$$\left\{\begin{array}{c}{\displaystyle {\mathfrak{D}}_{0^{+}}^{1.4}u\left(z\right)=u\left(z\right)+{\mathfrak{D}}_{0^{+}}^{0.3}u\left(z\right)+\widehat{\phi }\left(z\right),z\in (0,1),}\hfill \\ {\displaystyle u\left(0\right)=0,u\left(1\right)=8{\int }_0^{\frac{1}{2}}u\left(s\right)ds+54{\int }_0^{\frac{1}{3}}u\left(s\right)ds,}\hfill \end{array}\right.$$

(32)

where

$${\displaystyle \widehat{\phi }\left(z\right)=\frac{2}{\Gamma \left(1.6\right)}{z}^{0.6}-\frac{2}{\Gamma \left(2.7\right)}{z}^{1.7}-z^2.}$$

In this problem, we have taken data $\varrho =1.4$, $\xi =1/2$, $\eta =1/3$, ${\sigma }_n=0.3$, $\mu =\nu =1$, $p=8$ and $q=54$. It is known that $\varrho -{\sigma }_n=1.1>1$. In addition, $k_1(z,u\left(z\right))=k_2(z,u\left(z\right))=u\left(z\right)$ for $z\in [0,1]$. By assuming $m\left(z\right)={\mathfrak{D}}_{0^{+}}^{0.3}u\left(z\right)$, the equivalent integral equation of the problem (32) is the following

$$\begin{array}{cc}\hfill m\left(z\right)& ={\mathfrak{I}}_{0^{+}}^{1.1}\left[{\mathfrak{I}}_{0^{+}}^{0.3}m\left(z\right)+m\left(z\right)+\widehat{\phi }\left(z\right)\right]+\frac{\Gamma \left(1.4\right)}{\Gamma \left(1.1\right)}(8{\int }_0^{\frac{1}{2}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds\hfill \\ \hfill & +54{\int }_0^{\frac{1}{3}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds-{\mathfrak{I}}_{0^{+}}^{1.4}\left[{\mathfrak{I}}_{0^{+}}^{0.3}m\left(z\right)+m\left(z\right)+\widehat{\phi }\left(z\right)\right]{|}_{z=1})z^{0.1}\hfill \\ \hfill & ={\mathfrak{I}}_{0^{+}}^{1.4}m\left(z\right)+{\mathfrak{I}}_{0^{+}}^{1.1}m\left(z\right)+{\mathfrak{I}}_{0^{+}}^{1.1}\widehat{\phi }\left(z\right)+z^{0.1}\frac{8\Gamma \left(1.4\right)}{\Gamma \left(1.1\right)}{\int }_0^{\frac{1}{2}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds\hfill \\ \hfill & +z^{0.1}\frac{54\Gamma \left(1.4\right)}{\Gamma \left(1.1\right)}{\int }_0^{\frac{1}{3}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds-\frac{\Gamma \left(1.4\right)z^{0.1}}{\Gamma \left(1.1\right)}\left({\mathfrak{I}}_{0^{+}}^{1.7}m\left(z\right){|}_{z=1}\right)\hfill \\ \hfill & -\frac{\Gamma \left(1.4\right)z^{0.1}}{\Gamma \left(1.1\right)}\left({\mathfrak{I}}_{0^{+}}^{1.4}m\left(z\right){|}_{z=1}\right)-\frac{\Gamma \left(1.4\right)z^{0.1}}{\Gamma \left(1.1\right)}\left({\mathfrak{I}}_{0^{+}}^{1.4}\widehat{\phi }\left(z\right){|}_{z=1}\right).\hfill \end{array}$$

(33)

Thus, we decompose the right-hand side of (33) as

$$m\left(z\right)=\tilde{\mathbf{L}}\left(m\left(z\right)\right)+\tilde{\mathbf{N}}\left(m\left(z\right)\right)+\zeta \left(z\right),$$

where

$$\begin{array}{cc}\hfill \tilde{\mathbf{L}}\left(m\left(z\right)\right)& ={\mathfrak{I}}_{0^{+}}^{1.4}m\left(z\right)-{\mathfrak{I}}_{0^{+}}^{1.1}m\left(z\right),\hfill \end{array}$$

$$\begin{array}{cc}\hfill \tilde{\mathbf{N}}\left(m\left(z\right)\right)& =\frac{8\Gamma \left(1.4\right)z^{0.1}}{\Gamma \left(1.1\right)}{\int }_0^{\frac{1}{2}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds+\frac{54\Gamma \left(1.4\right)z^{0.1}}{\Gamma \left(1.1\right)}{\int }_0^{\frac{1}{3}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds\hfill \\ \hfill & -\frac{\Gamma \left(1.4\right)z^{0.1}}{\Gamma \left(1.1\right)}\left({\mathfrak{I}}_{0^{+}}^{1.7}m\left(z\right){|}_{z=1}\right)-\frac{\Gamma \left(1.4\right)z^{0.1}}{\Gamma \left(1.1\right)}\left({\mathfrak{I}}_{0^{+}}^{1.4}m\left(z\right){|}_{z=1}\right),\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill \zeta \left(z\right)& ={\mathfrak{I}}_{0^{+}}^{1.1}\widehat{\phi }\left(z\right)-\frac{\Gamma \left(1.4\right)z^{0.1}}{\Gamma \left(1.1\right)}\left({\mathfrak{I}}_{0^{+}}^{1.4}\widehat{\phi }\left(z\right){|}_{z=1}\right).\hfill \end{array}$$

Then, the sequence of approximate solutions of (32) and (33) are obtained by means of algorithms of the DGJIM and ADM methods as follows:

$••$Approximate solutions via DGJIM method for $\varrho =1.4$:

By using the suggested algorithm known as DGJIM numerical method in (24), we get

$$\begin{array}{cc}\hfill m_0\left(z\right)& =1.2948z^{1.7}-0.4262z^{2.8}-0.2936z^{3.1}-0.4748z^{0.1},\hfill \\ \hfill m_1\left(z\right)& =0.2936z^{3.1}-0.1228z^{4.2}-0.0382z^{4.5}-0.2398z^{1.5}+0.4261z^{2.8}-0.0968z^{3.9},\hfill \\ \hfill m_2\left(z\right)& =0.172z^{3.5}-0.0165z^{5.6}-0.0033z^{5.9}-0.0852z^{2.9}+0.1228z^{4.2}-0.0297z^{5.3}\hfill \\ \hfill & -0.0430z^{2.6}-0.0315z^{1.6}+0.0968z^{3.9}-0.0167z^5-0.1683z^{2.3}*0.4086z^{1.2}+4.4196z^{0.1}.\hfill \end{array}$$

Therefore,

$$\begin{array}{cc}\hfill w_0\left(z\right)& =1.2948z^{1.7}-0.4262z^{2.8}-0.2936z^{3.1}-0.4748z^{0.1},\hfill \\ \hfill w_1\left(z\right)& =1.2948z^{1.7}-0.1228z^{4.2}-0.0382z^{4.5}-0.3398z^{1.5}\hfill \\ \hfill & -0.0968z^{3.9}-0.4100z^{1.2}-4.4292z^{0.1},\hfill \\ \hfill w_2\left(z\right)& =1.2948z^{1.7}-0.0165z^{5.6}-0.0033z^{5.9}-0.0852z^{2.9}\hfill \\ \hfill & -0.0297z^{5.3}-0.0430z^{2.6}-0.0315z^{1.6}-0.0167z^5\hfill \\ \hfill & -0.0683z^{2.3}-0.0014z^{1.2}-0.0096z^{0.1}-0.0382z^{4.5}+0.2398z^{1.5},\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill u_0\left(z\right)& =z^2-0.2937z^{3.1}-0.1973z^{3.4}-0.5091z^{0.4},\hfill \\ \hfill u_1\left(z\right)& =z^2-0.0764z^{4.5}-0.0234z^{4.8}-0.2694z^{1.8}\hfill \\ \hfill & -0.0614z^{4.2}-0.3398z^{1.5}-4.7491z^{0.4},\hfill \\ \hfill u_2\left(z\right)& =z^2-0.0095z^{5.9}-0.0019z^{6.9}-0.0582z^{3.2}-0.0174z^{5.6}\hfill \\ \hfill & -0.0302z^{2.9}-0.0246z^{1.9}-0.0099z^{5.3}-0.0493z^{2.6}-0.0012z^{1.5}\hfill \\ \hfill & -0.0103z^{0.4}-0.0234z^{4.8}+0.1901z^{1.8}.\hfill \end{array}$$

$••$Approximate solutions via ADM method for $\varrho =1.4:$

By using the suggested algorithm known as ADM numerical method in (29), we get

$$\begin{array}{cc}\hfill m_0\left(z\right)& =1.2948z^{1.7}-0.4262z^{2.8}-0.2936z^{3.1}-0.4748z^{0.1},\hfill \\ \hfill m_1\left(z\right)& =0.2936z^{3.1}-0.1228z^{4.2}-0.0382z^{4.5}-0.3398z^{1.5}+0.4261z^{2.8}-0.0968z^{3.9}-0.4100z^{1.2}-3.9544z^{0.1},\hfill \\ \hfill m_2\left(z\right)& =0.1720z^{3.5}-0.0165z^{5.6}-0.0033z^{5.9}-0.0852z^{2.9}+0.1228z^{4.2}-0.0297z^{5.3}-0.2430z^{2.6}\hfill \\ \hfill & -0.0015z^{1.6}+0.0968z^{3.9}-0.0167z^5-0.1683z^{2.3}+0.4050z^{1.2}.\hfill \end{array}$$

Therefore,

$$\begin{array}{cc}\hfill w_0\left(z\right)& =1.2948z^{1.7}-0.4262z^{2.8}-0.2936z^{3.1}-0.4748z^{0.1},\hfill \\ \hfill w_1\left(z\right)& =1.2948z^{1.7}+4.4110z^{0.1}-0.1228z^{4.2}-0.0382z^{4.5}\hfill \\ \hfill & +0.3398z^{1.5}-0.0968z^{3.9}-0.4100z^{1.2},\hfill \\ \hfill w_2\left(z\right)& =1.2948z^{1.7}+0.1720z^{3.5}-0.0166z^{5.6}-0.0033z^{5.9}-0.0852z^{2.9}\hfill \\ \hfill & -0.0297z^{5.3}-0.2430z^{2.6}-0.0015z^{1.6}-0.0167z^5-0.1683z^{2.3}\hfill \\ \hfill & -0.0044z^{1.2}-0.0282z^{0.1}-0.0382z^{4.5}+0.3398z^{1.5},\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill u_0\left(z\right)& =z^2-0.0716z^{4.1}-0.1973z^{3.4}-0.5091z^{0.4},\hfill \\ \hfill u_1\left(z\right)& =z^2+4.7296z^{0.4}-0.0764z^{4.5}-0.0234z^{4.8}\hfill \\ \hfill & +0.2694z^{1.8}-0.0614z^{4.2}-0.3398z^{1.5},\hfill \\ \hfill u_2\left(z\right)& =z^2+0.1122z^{3.8}-0.0095z^{5.9}-0.0019z^{6.2}-0.0582z^{3.2}\hfill \\ \hfill & -0.0174z^{5.6}-0.1704z^{2.9}-0.0012z^{1.9}-0.0099z^{5.3}-0.1215z^{2.6}\hfill \\ \hfill & -0.0036z^{1.5}-0.0302z^{0.4}-0.0234z^{4.8}+0.2694z^{1.8}.\hfill \end{array}$$

In this case, the graphs of the three-term approximate solutions obtained by the DGJIM and ADM algorithms for the suggested RLFBVP (32) and the integral Equation (33) are plotted in Figure 1.

Note that, in view of Theorem 1, we prove that $u\left(z\right)$ is the solution of RLFBVP (1) if and only if $m\left(z\right)={\mathfrak{D}}_{0^{+}}^{{\sigma }_n}u\left(z\right)$ is the solution of the integral Equation (2). Now, in the case $\varrho =1.4$, since the exact solution of RLFBVP is given by $u\left(z\right)=z^2$, the corresponding exact solution of the equivalent integral equation is

$$m\left(z\right)={\mathfrak{D}}_{0^{+}}^{{\sigma }_n}{z}^2={\mathfrak{D}}_{0^{+}}^{0.3}{z}^2=\frac{2}{\Gamma \left(2.7\right)}{z}^{1.7}=1.2948z^{1.7}.$$

$•\mathbf{Case}\left(\mathbf{II}\right):\varrho =1.7$

In the next case, we consider the same problem for $\varrho =1.7$. In fact, at this time, we consider the following RLFBVP

$$\left\{\begin{array}{c}{\displaystyle {\mathfrak{D}}_{0^{+}}^{1.7}u\left(z\right)=u\left(z\right)+{\mathfrak{D}}_{0^{+}}^{0.3}u\left(z\right)+\widehat{\phi }\left(z\right),z\in (0,1),}\hfill \\ u\left(0\right)=0,\hfill \\ {\displaystyle u\left(1\right)=8{\int }_0^{\frac{1}{2}}u\left(s\right)ds+54{\int }_0^{\frac{1}{3}}u\left(s\right)ds,}\hfill \end{array}\right.$$

(34)

where

$${\displaystyle \widehat{\phi }\left(z\right)=\frac{2}{\Gamma \left(1.3\right)}{z}^{0.3}-\frac{2}{\Gamma \left(2.7\right)}{z}^{1.7}-z^2,}$$

such that we consider parameters $\varrho =1.7$, $\xi =1/2$, $\eta =1/3$, ${\sigma }_n=0.3$, $\mu =\nu =1$, $p=8$, and $q=54$. Obviously, $\varrho -{\sigma }_n=1.4>1$. In addition, $k_1(z,u\left(z\right))=k_2(z,u\left(z\right))=u\left(z\right)$ for $z\in [0,1]$. By assuming $m\left(z\right)={\mathfrak{D}}_{0^{+}}^{0.3}u\left(z\right)$, the equivalent integral equation of the problem (34) is given by

$$\begin{array}{cc}\hfill m\left(z\right)& ={\mathfrak{I}}_{0^{+}}^{1.4}\left[{\mathfrak{I}}_{0^{+}}^{0.3}m\left(z\right)+m\left(z\right)+\widehat{\phi }\left(z\right)\right]+\frac{\Gamma \left(1.7\right)}{\Gamma \left(1.4\right)}(8{\int }_0^{\frac{1}{2}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds\hfill \\ \hfill & +54{\int }_0^{\frac{1}{3}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds-{\mathfrak{I}}_{0^{+}}^{1.7}\left[{\mathfrak{I}}_{0^{+}}^{0.3}m\left(z\right)+m\left(z\right)+\widehat{\phi }\left(z\right)\right]{|}_{z=1})z^{0.4}\hfill \\ \hfill & ={\mathfrak{I}}_{0^{+}}^{1.7}m\left(z\right)+{\mathfrak{I}}_{0^{+}}^{1.4}m\left(z\right)+{\mathfrak{I}}_{0^{+}}^{1.4}\widehat{\phi }\left(z\right)+z^{0.4}\frac{8\Gamma \left(1.7\right)}{\Gamma \left(1.4\right)}{\int }_0^{\frac{1}{2}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds\hfill \\ \hfill & +z^{0.4}\frac{54\Gamma \left(1.7\right)}{\Gamma \left(1.4\right)}{\int }_0^{\frac{1}{3}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds-\frac{\Gamma \left(1.7\right)z^{0.4}}{\Gamma \left(1.4\right)}\left({\mathfrak{I}}_{0^{+}}^{2}m\left(z\right){|}_{z=1}\right)\hfill \\ \hfill & -\frac{\Gamma \left(1.7\right)z^{0.4}}{\Gamma \left(1.4\right)}\left({\mathfrak{I}}_{0^{+}}^{1.7}m\left(z\right){|}_{z=1}\right)-\frac{\Gamma \left(1.7\right)z^{0.4}}{\Gamma \left(1.4\right)}\left({\mathfrak{I}}_{0^{+}}^{1.7}\widehat{\phi }\left(z\right){|}_{z=1}\right).\hfill \end{array}$$

(35)

Then, we decompose the right-hand side of (35) as

$$m\left(z\right)=\tilde{\mathbf{L}}\left(m\left(z\right)\right)+\tilde{\mathbf{N}}\left(m\left(z\right)\right)+\zeta \left(z\right),$$

where

$$\begin{array}{cc}\hfill \tilde{\mathbf{L}}\left(m\left(z\right)\right)& ={\mathfrak{I}}_{0^{+}}^{1.7}m\left(z\right)+{\mathfrak{I}}_{0^{+}}^{1.4}m\left(z\right),\hfill \\ \hfill \tilde{\mathbf{N}}\left(m\left(z\right)\right)& =\frac{8\Gamma \left(1.7\right)z^{0.4}}{\Gamma \left(1.4\right)}{\int }_0^{\frac{1}{2}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds+\frac{54\Gamma \left(1.7\right)z^{0.4}}{\Gamma \left(1.4\right)}{\int }_0^{\frac{1}{3}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds\hfill \\ \hfill & -\frac{\Gamma \left(1.7\right)z^{0.4}}{\Gamma \left(1.4\right)}\left({\mathfrak{I}}_{0^{+}}^{2}m\left(z\right){|}_{z=1}\right)-\frac{\Gamma \left(1.7\right)z^{0.4}}{\Gamma \left(1.4\right)}\left({\mathfrak{I}}_{0^{+}}^{1.7}m\left(z\right){|}_{z=1}\right),\hfill \\ \hfill \zeta \left(z\right)& ={\mathfrak{I}}_{0^{+}}^{1.4}\widehat{\phi }\left(z\right)-\frac{\Gamma \left(1.7\right)z^{0.4}}{\Gamma \left(1.4\right)}\left({\mathfrak{I}}_{0^{+}}^{1.7}\widehat{\phi }\left(z\right){|}_{z=1}\right).\hfill \end{array}$$

Then, the sequence of approximate solutions of (34) and (35) are obtained by means of two DGJIM and ADM methods as follows:

$••$Approximate solutions via DGJIM method for $\varrho =1.7$:

$$\begin{array}{cc}\hfill w_0\left(z\right)& =1.2948z^{1.7}-0.3186z^{3.1}-0.1973z^{3.4}-0.6893z^{0.4},\hfill \\ \hfill w_1\left(z\right)& =1.2948z^{1.7}+0.0169z^{3.4}-2.7346z^{0.4}-0.0487z^{4.8}\hfill \\ \hfill & -0.1040z^{5.1}-0.2783z^{2.1}-0.1886z^{4.5}-0.3839z^{1.8},\hfill \\ \hfill w_2\left(z\right)& =1.2948z^{1.7}+0.0169z^{3.4}-0.0115z^{0.4}-0.0234z^{4.8}\hfill \\ \hfill & -0.0888z^{5.1}-0.0019z^{2.1}-0.1471z^{4.5}+0.2654z^{1.8}\hfill \\ \hfill & -0.0033z^{6.5}-0.1078z^{3.5}+0.0044z^{5.8}-0.0165z^{5.9},\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill u_0\left(z\right)& =z^2-0.2141z^{3.4}-0.1296z^{3.7}-0.6731z^{0.7},\hfill \\ \hfill u_1\left(z\right)& =z^2+0.0111z^{3.7}-2.6703z^{0.7}-0.1493z^{5.1}\hfill \\ \hfill & -0.0615z^{5.4}-0.2052z^{2.4}-0.0982z^{4.8}-0.2921z^{2.1},\hfill \\ \hfill u_2\left(z\right)& =z^2+0.0111z^{3.7}-0.0112z^{0.7}-0.0141z^{5.1}-0.0525z^{5.4}\hfill \\ \hfill & -0.0014z^{2.4}-0.0002z^{7.1}-0.0449z^{4.1}-0.0075z^{6.5}-0.0703z^{3.8}\hfill \\ \hfill & +0.0025z^{6.1}-0.0018z^{6.8}-0.0094z^{6.2}-0.0899z^{4.8}+0.2025z^{2.1}.\hfill \end{array}$$

$••$Approximate solutions via ADM method for $\varrho =1.7$:

$$\begin{array}{cc}\hfill w_0\left(z\right)& =1.2948z^{1.7}-0.3186z^{3.1}-0.1973z^{3.4}-0.6893z^{0.4},\hfill \\ \hfill w_1\left(z\right)& =1.2948z^{1.7}+0.0169z^{3.4}-0.0346z^{0.4}-0.0487z^{4.8}\hfill \\ \hfill & -0.1040z^{5.1}-0.2783z^{2.1}-0.1886z^{4.5}-0.3839z^{1.8},\hfill \\ \hfill w_2\left(z\right)& =1.2948z^{1.7}+0.0169z^{3.4}-0.0346z^{0.4}-0.0234z^{4.8}+0.0012z^{5.1}\hfill \\ \hfill & +1.7119z^{2.1}-0.1471z^{4.5}-1.3654z^{1.8}-0.0033z^{6.5}-0.0005z^{6.8}\hfill \\ \hfill & -0.0703z^{3.8}-0.0134z^{6.2}-0.1078z^{3.5}+0.0044z^{5.8}-0.0165z^{5.9},\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill u_0\left(z\right)& =z^2-0.2141z^{3.4}-0.1296z^{3.7}-0.6731z^{0.7},\hfill \\ \hfill u_1\left(z\right)& =z^2+0.1055z^{3.7}-0.0338z^{0.7}-0.0111z^{3.7}-0.0293z^{5.1}\hfill \\ \hfill & -0.0083z^{5.4}-0.2052z^{2.4}-0.1153z^{4.8}-0.2921z^{2.1},\hfill \\ \hfill u_2\left(z\right)& =z^2+0.0111z^{3.7}-0.0338z^{0.7}+0.0111z^{3.7}-0.0141z^{5.1}\hfill \\ \hfill & +0.0007z^{5.4}+1.2619z^{2.4}-0.0899z^{4.8}-1.0416z^{2.1}-0.0018z^{6.8}\hfill \\ \hfill & -0.0002z^{7.1}-0.0449z^{4.1}-0.0075z^{6.5}-0.0703z^{3.8}+0.0025z^{6.1}-0.0094z^{6.2}.\hfill \end{array}$$

In consequence, the graphs of the three-term approximate solutions obtained by the DGJIM and ADM algorithm for the suggested RLFBVP (34) and the integral Equation (35) are plotted in Figure 2.

$•\mathbf{Case}\left(\mathbf{III}\right):\varrho =1.9$

Finally, we consider the first problem for $\varrho =1.9$ as the third case. Consider the following RLFBVP

$$\left\{\begin{array}{c}{\displaystyle {\mathfrak{D}}_{0^{+}}^{1.9}u\left(z\right)=u\left(z\right)+{\mathfrak{D}}_{0^{+}}^{0.3}u\left(z\right)+\widehat{\phi }\left(z\right),z\in (0,1),}\hfill \\ u\left(0\right)=0,\hfill \\ {\displaystyle u\left(1\right)=8{\int }_0^{\frac{1}{2}}u\left(s\right)ds+54{\int }_0^{\frac{1}{3}}u\left(s\right)ds,}\hfill \end{array}\right.$$

(36)

where

$${\displaystyle \widehat{\phi }\left(z\right)=\frac{2}{\Gamma \left(1.1\right)}{z}^{0.1}-\frac{2}{\Gamma \left(2.7\right)}{z}^{1.7}-z^2.}$$

The parameters $\varrho =1.9$, $\xi =1/2$, $\eta =1/3$, ${\sigma }_n=0.3$, $\mu =\nu =1$, $p=8$, and $q=54$ are assumed here. Evidently, $\varrho -{\sigma }_n=1.6>1$. In addition, $k_1(z,u\left(z\right))=k_2(z,u\left(z\right))=u\left(z\right)$ for $z\in [0,1]$. By assuming $m\left(z\right)={\mathfrak{D}}_{0^{+}}^{0.3}u\left(z\right)$, the equivalent integral equation of the problem (36) is given in the following form

$$\begin{array}{cc}\hfill m\left(z\right)& ={\mathfrak{I}}_{0^{+}}^{1.6}\left[{\mathfrak{I}}_{0^{+}}^{0.3}m\left(z\right)+m\left(z\right)+\widehat{\phi }\left(z\right)\right]+\frac{\Gamma \left(1.9\right)}{\Gamma \left(1.6\right)}(8{\int }_0^{\frac{1}{2}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds\hfill \\ \hfill & +54{\int }_0^{\frac{1}{3}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds-{\mathfrak{I}}_{0^{+}}^{1.9}\left[{\mathfrak{I}}_{0^{+}}^{0.3}m\left(z\right)+m\left(z\right)+\widehat{\phi }\left(z\right)\right]{|}_{z=1})z^{0.6}\hfill \\ \hfill & ={\mathfrak{I}}_{0^{+}}^{1.9}m\left(z\right)+{\mathfrak{I}}_{0^{+}}^{1.6}m\left(z\right)+{\mathfrak{I}}_{0^{+}}^{1.6}\widehat{\phi }\left(z\right)+z^{0.6}\frac{8\Gamma \left(1.9\right)}{\Gamma \left(1.6\right)}{\int }_0^{\frac{1}{2}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds\hfill \\ \hfill & +z^{0.6}\frac{54\Gamma \left(1.9\right)}{\Gamma \left(1.6\right)}{\int }_0^{\frac{1}{3}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds-\frac{\Gamma \left(1.9\right)z^{0.6}}{\Gamma \left(1.6\right)}\left({\mathfrak{I}}_{0^{+}}^{2.2}m\left(z\right){|}_{z=1}\right)\hfill \\ \hfill & -\frac{\Gamma \left(1.9\right)z^{0.6}}{\Gamma \left(1.6\right)}\left({\mathfrak{I}}_{0^{+}}^{1.9}m\left(z\right){|}_{z=1}\right)-\frac{\Gamma \left(1.9\right)z^{0.6}}{\Gamma \left(1.6\right)}\left({\mathfrak{I}}_{0^{+}}^{1.9}\widehat{\phi }\left(z\right){|}_{z=1}\right).\hfill \end{array}$$

(37)

By decomposing the right-hand side of (37), we get

$$m\left(z\right)=\tilde{\mathbf{L}}\left(m\left(z\right)\right)+\tilde{\mathbf{N}}\left(m\left(z\right)\right)+\zeta \left(z\right),$$

where

$$\begin{array}{cc}\hfill \tilde{\mathbf{L}}\left(m\left(z\right)\right)& ={\mathfrak{I}}_{0^{+}}^{1.9}m\left(z\right)+{\mathfrak{I}}_{0^{+}}^{1.6}m\left(z\right),\hfill \\ \hfill \tilde{\mathbf{N}}\left(m\left(z\right)\right)& =\frac{8\Gamma \left(1.9\right)z^{0.6}}{\Gamma \left(1.6\right)}{\int }_0^{\frac{1}{2}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds+\frac{54\Gamma \left(1.9\right)z^{0.6}}{\Gamma \left(1.6\right)}{\int }_0^{\frac{1}{3}}{\mathfrak{I}}_{0^{+}}^{0.3}m\left(s\right)ds\hfill \\ \hfill & -\frac{\Gamma \left(1.9\right)z^{0.6}}{\Gamma \left(1.6\right)}\left({\mathfrak{I}}_{0^{+}}^{2.2}m\left(z\right){|}_{z=1}\right)-\frac{\Gamma \left(1.9\right)z^{0.6}}{\Gamma \left(1.6\right)}\left({\mathfrak{I}}_{0^{+}}^{1.9}m\left(z\right){|}_{z=1}\right),\hfill \\ \hfill \zeta \left(z\right)& ={\mathfrak{I}}_{0^{+}}^{1.6}\widehat{\phi }\left(z\right)-\frac{\Gamma \left(1.9\right)z^{0.6}}{\Gamma \left(1.6\right)}\left({\mathfrak{I}}_{0^{+}}^{1.9}\widehat{\phi }\left(z\right){|}_{z=1}\right).\hfill \end{array}$$

Then, the sequence of approximate solutions are obtained by means of two DGJIM and ADM methods illustrated as:

$••$Approximate solutions via DGJIM method for $\varrho =1.9$:

$$\begin{array}{cc}\hfill w_0\left(z\right)& =1.2948z^{1.7}-0.2259z^{3.3}-0.1495z^{3.6}-0.8114z^{0.6},\hfill \\ \hfill w_1\left(z\right)& =1.2948z^{1.7}-1.8726z^{0.6}-0.0236z^{5.2}-0.0069z^{5.5}\hfill \\ \hfill & -0.2182z^{2.5}-0.0198z^{4.9}-0.2991z^{2.2},\hfill \\ \hfill w_2\left(z\right)& =1.2948z^{1.7}+0.0427z^{0.6}+0.0017z^{5.2}-0.5008z^{2.5}\hfill \\ \hfill & +0.4866z^{2.2}-0.0009z^{7.1}-0.0001z^{7.4}-0.0163z^{4.4}\hfill \\ \hfill & -0.0017z^{6.8}-0.0520z^{4.1}-0.0011z^{6.5}-0.0406z^{3.8},\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill u_0\left(z\right)& =z^2-0.1495z^{3.6}-0.0968z^{3.9}-0.7538z^{0.9},\hfill \\ \hfill u_1\left(z\right)& =z^2-1.7397z^{0.9}-0.0139z^{5.5}-0.0040z^{5.8}\hfill \\ \hfill & -0.1545z^{2.8}-0.0118z^{5.2}-0.2182z^{2.5},\hfill \\ \hfill u_2\left(z\right)& =z^2+0.0397z^{0.9}+0.0010z^{5.5}-0.3546z^{2.8}\hfill \\ \hfill & +0.3549z^{25}-0.0004z^{7.4}-0.0005z^{7.7}-0.0118z^{4.7}\hfill \\ \hfill & -0.0009z^{7.1}-0.0326z^{4.4}-0.0006z^{6.8}-0.0025z^{4.1}.\hfill \end{array}$$

$••$Approximate solutions via ADM method for $\varrho =1.9$:

$$\begin{array}{cc}\hfill w_0\left(z\right)& =1.2948z^{1.7}-0.2259z^{3.3}-0.1495z^{3.6}-0.8114z^{0.6},\hfill \\ \hfill w_1\left(z\right)& =1.2948z^{1.7}-1.8626z^{0.6}-0.0236z^{5.2}-0.0069z^{5.5}\hfill \\ \hfill & -0.2182z^{2.5}-0.0198z^{4.9}-0.2991z^{2.2},\hfill \\ \hfill w_2\left(z\right)& =1.2948z^{1.7}-0.0126z^{0.6}+0.0017z^{5.2}-0.5008z^{2.5}\hfill \\ \hfill & +0.4866z^{2.2}-0.0009z^{7.1}-0.0001z^{7.4}-0.0163z^{4.4}\hfill \\ \hfill & -0.0017z^{6.8}-0.0520z^{4.1}-0.0011z^{6.5}-0.0406z^{3.8},\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill u_0\left(z\right)& =z^2-0.1495z^{3.6}-0.0968z^{3.9}-0.7528z^{0.9},\hfill \\ \hfill u_1\left(z\right)& =z^2-1.7304z^{0.9}-0.0139z^{5.5}-0.0040z^{5.8}\hfill \\ \hfill & -0.1545z^{2.8}-0.0118z^{5.2}-0.2182z^{2.5},\hfill \\ \hfill u_2\left(z\right)& =z^2-0.0117z^{0.9}+0.0010z^{5.5}-0.3546z^{2.8}\hfill \\ \hfill & +0.3549z^{2.5}-0.0004z^{7.41}-0.0005z^{7.7}-0.0100z^{4.7}\hfill \\ \hfill & -0.0009z^{7.1}-0.0626z^{4.4}-0.0006z^{6.8}-0.0259z^{4.1}.\hfill \end{array}$$

In consequence, the graphs of the three-term approximate solutions obtained by the DGJIM and ADM algorithm for the suggested RLFBVP (36) and the integral Equation (37) are plotted in Figure 3.

6. Conclusions

In this paper, we study the existence of solutions for a multi-term multi-order RLFBVP with integral boundary conditions in the first step. Next, we apply two numerical methods (i.e., DGJIM and ADM algorithms) for solving the suggested multi-term fractional differential equation based on the decomposition technique. We show by an example that the approximate solutions obtained by these methods are in excellent agreement with the exact solutions. These give the solution as a series that quickly converges to the exact one if it exists. Therefore, this paper states that these two numerical methods can be utilized in many other multi-term FBVPs with different boundary value conditions by terms of some symmetric and asymmetric operators.