H-U-Type Stability and Numerical Solutions for a Nonlinear Model of the Coupled Systems of Navier BVPs via the Generalized Differential Transform Method

Abstract

This paper is devoted to generalizing the standard system of Navier boundary value problems to a fractional system of coupled sequential Navier boundary value problems by using terms of the Caputo derivatives. In other words, for the first time, we design a multi-term fractional coupled system of Navier equations under the fractional boundary conditions. The existence theory is studied regarding solutions of the given coupled sequential Navier boundary problems via the Krasnoselskii’s fixed-point theorem on two nonlinear operators. Moreover, the Banach contraction principle is applied to investigate the uniqueness of solution. We then focus on the Hyers–Ulam-type stability of its solution. Furthermore, the approximate solutions of the proposed coupled fractional sequential Navier system are obtained via the generalized differential transform method. Lastly, the results of this research are supported by giving simulated examples.

1. Introduction

Fractional differential equations (FDEs) are considered an important area of research in the direction of the applications of fractional calculus. Regarding applications of FDEs, one can observe some fields including aerodynamics, biochemistry, electro-chemistry, bioengineering, physics, viscoelasticity, mathematical biology, and so on (see [1,2]). In the last few years, a large number of studies regarding the existence theory for different FDEs have received much attentions from researchers and some examples include [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]. Since we can model some of applied phenomena in the framework of the fractional coupled systems, a large number of researchers have conducted many research studies on the existence of solution for such a type of systems (for instances, see [19,20,21,22,23]).

Stability analysis along with numerical techniques are the most important components of research in this regard. Usually, obtaining the exact solutions of a non-linear boundary problem in the fractional settings is time-consuming wor,k and it is a task full of challenges. Therefore, these items motivated mathematicians to find the best approximate solutions for existing boundary value problems (BVPs). To perform this, various procedures were introduced, including decomposition methods [24,25], HATM [26], q-HATM [27], and integral transforms [28].

The generalized differential transform method (GDT-method) is considered as one of the most accurate and strongest techniques for finding approximate solutions of nonlinear or linear FDEs. This transform can be observed in some papers such as [29,30] in which the authors have used it to analyze approximate solutions of the given BVPs. Notice that there is no specified method to study nonlinear FDEs for obtaining explicit solutions. Therefore, a reliable method is required to find approximate solutions in the framework of finite/infinite series regarding the given FBVPs.

In 2018, Shah et al. [31] established the existence results and obtained approximate solutions with the help of GDT-method for the system of coupled FDEs with movable integral conditions:

$$\left\{\begin{array}{c}{}^c{\mathfrak{D}}^p\mu \left(\mathfrak{t}\right)=y_1(\mathfrak{t},\nu \left(\mathfrak{t}\right)),(\forall \mathfrak{t}\in \mathbb{I}),\hfill \\ {}^c{\mathfrak{D}}^q\nu \left(\mathfrak{t}\right)=y_2(\mathfrak{t},\mu \left(\mathfrak{t}\right)),(\forall \mathfrak{t}\in \mathbb{I}),\hfill \\ \mu \left(0\right)=0,\nu \left(0\right)=0,\hfill \\ {\displaystyle \mu \left(1\right)={\int }_0^{c_1}\mu \left(r\right)\mathrm{d}r,\nu \left(1\right)={\int }_0^{c_2}\nu \left(r\right)\mathrm{d}r,}\hfill \end{array}\right.$$

by assuming $1<p,q\le 2$, $y_1,y_2\in C(\mathbb{I}\times [0,\infty ))$, $\mathbb{I}:=[0,1]$, and $0<c_1,c_2<1$. In 2020, Alrabaiah, Ahmad, Shah, and Rahman [32] conducted qualitative research regarding H-U-type stability of solutions to a system of coupled nonlinear integral delay pantograph BVPs:

$$\left\{\begin{array}{c}{\mathfrak{D}}_0^{{\ell }_1}\mu \left(\mathfrak{t}\right)+{\chi }_1(\mathfrak{t},\mu \left(a\mathfrak{t}\right),\nu \left(\mathfrak{t}\right),{\mathfrak{D}}_0^{q_1}\nu \left(\mathfrak{t}\right))=0,(\forall \mathfrak{t}\in \mathbb{I}),\hfill \\ {\mathfrak{D}}_0^{{\ell }_2}\nu \left(\mathfrak{t}\right)+{\chi }_2(\mathfrak{t},\mu \left(\mathfrak{t}\right),\nu \left(a\mathfrak{t}\right),{\mathfrak{D}}_0^{q_2}\mu \left(\mathfrak{t}\right))=0,\hfill \\ {\displaystyle \mu \left(0\right)=0,\mu \left(1\right)={\int }_0^{1}y\left(r\right)\mu \left(r\right)\mathrm{d}r,}\hfill \\ {\displaystyle \nu \left(0\right)=0,\nu \left(1\right)={\int }_0^{1}y\left(r\right)\nu \left(r\right)\mathrm{d}r,}\hfill \end{array}\right.$$

with $1<{\ell }_1,{\ell }_2\le 2$. $\mathbb{I}:=[0,1]$, $0<q_1,q_2,a<1$, and ${\chi }_1,{\chi }_2:\mathbb{I}\times {\mathbb{R}}^3\to \mathbb{R}$ are nonlinear, and $y:(0,1)\to [0,\infty )$ is bounded.

In this paper, inspired by [33] and the above papers, we focus on the intention in which some qualitative aspects of possible solutions for a system of the coupled fractional sequential Navier model are investigated. In more precise words, we consider the following structure of a coupled sequential Navier FBVPs as follows:

$$\left\{\begin{array}{c}{}^C{\mathfrak{D}}^{{\ell }_1}\left({}^C{\mathfrak{D}}^{{\ell }_2}\mu \right)\left(\mathfrak{t}\right)=K\left(\mathfrak{t},\nu \left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(\mathfrak{t}\right)\right),(\mathfrak{t}\in \mathcal{O}:=[0,1]),\hfill \\ {}^C{\mathfrak{D}}^{{\ell }_1^{*}}\left({}^C{\mathfrak{D}}^{{\ell }_2^{*}}\nu \right)\left(\mathfrak{t}\right)=M\left(\mathfrak{t},\mu \left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(\mathfrak{t}\right)\right),(\mathfrak{t}\in \mathcal{O}:=[0,1]),\hfill \\ \gamma \mu \left(0\right)=\delta \mu \left(1\right)=\alpha {}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(0\right)=\beta {}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(1\right)=0,\hfill \\ {\gamma }^{*}\nu \left(0\right)={\delta }^{*}\nu \left(1\right)={\alpha }^{*}{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\nu \left(0\right)={\beta }^{*}{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\nu \left(1\right)=0,\hfill \end{array}\right.$$

(1)

where ${\ell }_1,{\ell }_1^{*}\in (1,2]$, ${\ell }_2,{\ell }_2^{*}\in (1,2]$ and $\gamma ,\delta ,\alpha ,\beta ,{\gamma }^{*},{\delta }^{*},{\alpha }^{*},{\beta }^{*}\in {\mathbb{R}}^{+}$. Moreover, the operator ${}^C{\mathfrak{D}}^{(·)}$ denotes the Caputo fractional derivative of different orders. Furthermore, two continuous single-valued functions $K,M:\mathcal{O}\times {\mathbb{R}}^2\to \mathbb{R}$ are assumed to be arbitrary equipped with some needed properties, which are explained in the sequel.

Notice that the novelty of our paper is that the above suggested structure for the Navier problem is unique and novel, and one can consider it as a generalized fractional model of standard Navier problem in the context of Caputo operators. Indeed, by taking ${\ell }_1={\ell }_1^{*}={\ell }_2={\ell }_2^{*}=2$ and $\gamma ={\gamma }^{*}=\delta ={\delta }^{*}=\alpha ={\alpha }^{*}=\beta ={\beta }^{*}=1$, we obtain the standard system of Navier BVPs of the fourth orders formulated by the following.

$$\left\{\begin{array}{c}{\mu }^{\left(4\right)}\left(\mathfrak{t}\right)=K\left(\mathfrak{t},\nu \left(\mathfrak{t}\right),{\nu }^{″}\left(\mathfrak{t}\right)\right),(\mathfrak{t}\in \mathcal{O}:=[0,1]),\hfill \\ {\nu }^{\left(4\right)}\left(\mathfrak{t}\right)=M\left(\mathfrak{t},\mu \left(\mathfrak{t}\right),{\mu }^{″}\left(\mathfrak{t}\right)\right),(\mathfrak{t}\in \mathcal{O}:=[0,1]),\hfill \\ \mu \left(0\right)=\mu \left(1\right)={\mu }^{″}\left(0\right)={\mu }^{″}\left(1\right)=0,\hfill \\ \nu \left(0\right)=\nu \left(1\right)={\nu }^{″}\left(0\right)={\nu }^{″}\left(1\right)=0.\hfill \end{array}\right.$$

(2)

For deriving the needed conditions implying the uniqueness and existence of solutions to a given coupled system of sequential Navier FBVPs, fixed point theory is used in this paper. In this direction, the Krasnoselskii’s fixed point theorem and Banach’s contraction principle have the key role for proving the existence and uniqueness of solution. To study the H-U-type stability criterion, we first define the relevant notions for the given system (Equation (1)) and then prove the results with respect to this qualitative property. Along with these, in the sequel, the approximate solutions of Equation (1) are obtained by applying numerical algorithms arising in the GDT-method. All these approximate solutions are plotted in two illustrative examples. We emphasize that the main motivation is that we conduct our research based on numerical methods of differential transform type to search for approximate solutions of a new fractional model of Navier BVP, which is more applicable in fourth-order models arising in engineering, and this renders our theoretical findings as useful results. In other words, for the first time, we apply the generalized differential transform method for a real model of boundary value problems to investigate the qualitative behaviors of the given system. Until now, the existing limited works regarding the Navier problem studied the existence results for integer-order systems, while we not only generalize the supposed system to a fractional structure but we also try to use an accurate algorithm for finding the approximate solution. These results show the novelty of our research. By the accuracy of this method, we can analyze other models in this direction.

The rest of the contents are as follows: Preliminaries are collected in the next section. The results in relation to the existence of solution and its uniqueness are presented in Section 3. In Section 4, the results regarding H-U-type stability criterion are proved. Section 5 introduces numerical algorithms of the GDT-method for a given sequential system (1). Different cases of the sequential Navier system are analyzed in two examples (with graphs) in Section 6. The last section is devoted to presenting conclusive remarks.

2. Preliminaries

As we will observe, two concepts of the fractional operators such as the Riemann–Liouville integral and the Caputo derivative play useful roles in this study. Therefore, we recall several properties of them here.

Definition 1

([2,34]). The Riemann–Liouville integral operator of the fractional order ℓ for a given function $\mu :[0,+\infty )\to \mathbb{R}$ is defined by the following:

$${}^R{\mathfrak{I}}^{\ell }\mu \left(\mathfrak{t}\right)=\frac{1}{\mathsf{\Gamma }(\ell )}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{\ell -1}\mu \left(r\right)\mathrm{d}r,\ell >0,$$

if it exists.

Definition 2

([2,34]). By taking $m=1+[\ell ]$, for $\mu \in C({\mathbb{R}}^{\ge 0},\mathbb{R})$, the ℓ-th-Riemann–Liouville derivation operator is defined by the following:

$${\displaystyle {}^R{\mathfrak{D}}^{\ell }\mu \left(\mathfrak{t}\right)={\left(\frac{\mathrm{d}}{\mathrm{d}\mathfrak{t}}\right)}^m\frac{1}{\mathsf{\Gamma }(m-\ell )}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{m-\ell -1}\mu \left(r\right)\mathrm{d}r,}$$

if it exists.

Definition 3

([2,34]). By taking $m=1+[\ell ]$, for $\mu \in A{C}^{\left(m\right)}({\mathbb{R}}^{\ge 0},\mathbb{R})$, the ℓ-th-Caputo derivation operator is defined by the following:

$${\displaystyle {}^C{\mathfrak{D}}^{\ell }\mu \left(\mathfrak{t}\right)=\frac{1}{\mathsf{\Gamma }(m-\ell )}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{m-\ell -1}{\mu }^{\left(m\right)}\left(r\right)\mathrm{d}r,}$$

if it exists.

Proposition 1

([1]). If $\ell \in (m-1,m)$, then for each $\mu \in C^{m-1}(0,\infty )$:

$${}^R{\mathfrak{I}}^{\ell }\left({}^C{\mathfrak{D}}^{\ell }\mu \right)\left(\mathfrak{t}\right)=\mu \left(\mathfrak{t}\right)+{\mathfrak{c}}_0+{\mathfrak{c}}_1\mathfrak{t}+{\mathfrak{c}}_2{\mathfrak{t}}^2+\cdots +{\mathfrak{c}}_{m-1}{\mathfrak{t}}^{m-1},$$

is valid for some ${\mathfrak{c}}_0,{\mathfrak{c}}_1,\cdots ,{\mathfrak{c}}_{m-1}\in \mathbb{R}$.

In the next section, we establish some results regarding the existence criteria.

3. Results for Existence

It is an evident notion that $\mathfrak{B}=\{\mu \left(\mathfrak{t}\right):\mu \left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(\mathfrak{t}\right)\in C_{{\mathbb{R}}^{\ge 0}}\left(\mathcal{O}\right)\}$ is a space of the Banach type subject to the norm ${\parallel \mu \parallel }_{\mathfrak{B}}={max}_{\mathfrak{t}\in \mathcal{O}}|\mu \left(\mathfrak{t}\right)|+{max}_{\mathfrak{t}\in \mathcal{O}}|{}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(\mathfrak{t}\right)|$. Consequently, $\mathfrak{B}\times \mathfrak{B}$ is a product Banach space that is equipped with norm ${\parallel (\mu ,\nu )\parallel }_{\mathfrak{B}\times \mathfrak{B}}={max\{\parallel \mu \parallel }_{\mathfrak{B}}{,\parallel \nu \parallel }_{\mathfrak{B}}\}$.

Proposition 2.

Let $T\in C_{\mathbb{R}}\left(\mathcal{O}\right)$, ${\ell }_1\in (1,2)$, ${\ell }_2\in (1,2)$, and $\gamma ,\delta ,\alpha ,\beta \in {\mathbb{R}}^{+}$. Then the solution of the nonlinear sequential fractional Navier BVP defined by the following:

$$\left\{\begin{array}{c}{}^C{\mathfrak{D}}^{{\ell }_1}\left({}^C{\mathfrak{D}}^{{\ell }_2}\mu \right)\left(\mathfrak{t}\right)=T\left(\mathfrak{t}\right),(\mathfrak{t}\in \mathcal{O}),\hfill \\ \gamma \mu \left(0\right)=\delta \mu \left(1\right)=\alpha {}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(0\right)=\beta {}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(1\right)=0,\hfill \end{array}\right.$$

(3)

is given by

$$\begin{array}{cc}\hfill \mu \left(\mathfrak{t}\right)& ={\int }_0^{\mathfrak{t}}\frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}T\left(r\right)\mathrm{d}r-\mathfrak{t}{\int }_0^1\frac{{(1-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}T\left(r\right)\mathrm{d}r\hfill \\ \hfill & +\frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}T\left(r\right)\mathrm{d}r.\hfill \end{array}$$

(4)

Proof. 

We consider the function $\mu$ as a solution for the nonlinear sequential Navier BVP (3). Then, we obtain the following: ${}^C{\mathfrak{D}}^{{\ell }_1}\left({}^C{\mathfrak{D}}^{{\ell }_2}\mu \right)\left(\mathfrak{t}\right)=T\left(\mathfrak{t}\right)$. By virtue of ${\ell }_1\in (1,2)$ and by taking the fractional integral in the Riemann-Liouville settings of order ${\ell }_1$ on both sides of the latter equation, we get

$${}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(\mathfrak{t}\right)={\int }_0^{\mathfrak{t}}\frac{{(\mathfrak{t}-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}T\left(r\right)\mathrm{d}r+{\mathfrak{c}}_0+{\mathfrak{c}}_1\mathfrak{t},$$

so that ${\mathfrak{c}}_0,{\mathfrak{c}}_1\in \mathbb{R}$ are some constants that we need to find. Now, by the third condition $\alpha {}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(0\right)=0$, we obtain ${\mathfrak{c}}_0=0$. Thus, the following is the case.

$${}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(\mathfrak{t}\right)={\int }_0^{\mathfrak{t}}\frac{{(\mathfrak{t}-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}T\left(r\right)\mathrm{d}r+{\mathfrak{c}}_1\mathfrak{t}.$$

(5)

On the other hand, by considering Equation (5) and the fourth condition $\beta {}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(1\right)=0$, we obtain the following.

$$\beta {\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}T\left(r\right)\mathrm{d}r+\beta {\mathfrak{c}}_1=0,$$

Thus, the following is the case.

$${\mathfrak{c}}_1=-{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}T\left(r\right)\mathrm{d}r.$$

(6)

In the sequel, in view of Equation (6), relation (5) becomes the following.

$${}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(\mathfrak{t}\right)={\int }_0^{\mathfrak{t}}\frac{{(\mathfrak{t}-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}T\left(r\right)\mathrm{d}r-\mathfrak{t}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}T\left(r\right)\mathrm{d}r.$$

(7)

Again, by virtue of ${\ell }_2\in (1,2)$ and by taking the fractional integral in the Riemann–Liouville settings of order ${\ell }_2$ on both sides of the latter equation, we obtain the following:

$$\mu \left(\mathfrak{t}\right)={\int }_0^{\mathfrak{t}}\frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}T\left(r\right)\mathrm{d}r-\frac{{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}T\left(r\right)\mathrm{d}r+{\mathfrak{c}}_0^{*}+{\mathfrak{c}}_1^{*}\mathfrak{t},$$

in which ${\mathfrak{c}}_0^{*}$ and ${\mathfrak{c}}_1^{*}\in \mathbb{R}$ are some constants that we have to find. Immediately, the first condition $\gamma \mu \left(\mathfrak{0}\right)=0$ gives us ${\mathfrak{c}}_0^{*}=0$. Consequently, we have the following.

$$\mu \left(\mathfrak{t}\right)={\int }_0^{\mathfrak{t}}\frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}T\left(r\right)\mathrm{d}r-\frac{{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}T\left(r\right)\mathrm{d}r+{\mathfrak{c}}_1^{*}\mathfrak{t}.$$

(8)

Lastly, the second condition $\delta \mu \left(1\right)=0$ implies that the following is the case.

$$\delta {\int }_0^1\frac{{(1-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}T\left(r\right)\mathrm{d}r-\frac{\delta }{\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}T\left(r\right)\mathrm{d}r+\delta {\mathfrak{c}}_1^{*}=0.$$

Consequently, the constant ${\mathfrak{c}}_1^{*}$ is obtained as follows.

$${\mathfrak{c}}_1^{*}=-{\int }_0^1\frac{{(1-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}T\left(r\right)\mathrm{d}r+\frac{1}{\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}T\left(r\right)\mathrm{d}r.$$

We insert ${\mathfrak{c}}_1^{*}$ into Equation (8) and obtain the following.

$$\begin{array}{cc}\hfill \mu \left(\mathfrak{t}\right)& ={\int }_0^{\mathfrak{t}}\frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}T\left(r\right)\mathrm{d}r-\mathfrak{t}{\int }_0^1\frac{{(1-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}T\left(r\right)\mathrm{d}r\hfill \\ \hfill & +\frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}T\left(r\right)\mathrm{d}r,\hfill \end{array}$$

This yields the result that $\mu$ satisfies the integral Equation (4), and the proof is completed. □

In view of above proposition, we here present an equivalent version of the coupled system of integral equations in relation to the given system of coupled sequential Navier BVPs (1) by the following proposition.

Proposition 3.

Let ${\ell }_1,{\ell }_1^{*}\in (1,2]$, ${\ell }_2,{\ell }_2^{*}\in (1,2]$ and $\gamma ,\delta ,\alpha ,\beta ,{\gamma }^{*},{\delta }^{*},{\alpha }^{*},{\beta }^{*}\in {\mathbb{R}}^{+}$ and $K,M\in C_{{\mathbb{R}}^{\ge 0}}$$(\mathcal{O}\times {\mathbb{R}}^{\ge 0}\times {\mathbb{R}}^{\ge 0})$. Then, an equivalent version of the coupled BVPs of fractional Navier differential equations is as follows:

$$\left\{\begin{array}{c}{}^C{\mathfrak{D}}^{{\ell }_1}\left({}^C{\mathfrak{D}}^{{\ell }_2}\mu \right)\left(\mathfrak{t}\right)=K\left(\mathfrak{t},\nu \left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(\mathfrak{t}\right)\right),(\mathfrak{t}\in \mathcal{O}:=[0,1]),\hfill \\ {}^C{\mathfrak{D}}^{{\ell }_1^{*}}\left({}^C{\mathfrak{D}}^{{\ell }_2^{*}}\nu \right)\left(\mathfrak{t}\right)=M\left(\mathfrak{t},\mu \left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(\mathfrak{t}\right)\right),(\mathfrak{t}\in \mathcal{O}:=[0,1]),\hfill \\ \gamma \mu \left(0\right)=\delta \mu \left(1\right)=\alpha {}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(0\right)=\beta {}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(1\right)=0,\hfill \\ {\gamma }^{*}\nu \left(0\right)={\delta }^{*}\nu \left(1\right)={\alpha }^{*}{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\nu \left(0\right)={\beta }^{*}{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\nu \left(1\right)=0,\hfill \end{array}\right.$$

(9)

This is given by the following coupled integral equations.

$$\left\{\begin{array}{c}{\displaystyle \mu \left(\mathfrak{t}\right)={\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)\mathrm{d}r-\mathfrak{t}{\int }_0^1{\displaystyle \frac{{(1-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)\mathrm{d}r}\hfill \\ {\displaystyle +{\displaystyle \frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+2)}}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)\mathrm{d}r,(\forall \mathfrak{t}\in \mathcal{O}),}\hfill \\ {\displaystyle \nu \left(\mathfrak{t}\right)={\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1^{*}+{\ell }_2^{*}-1}}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*})}}M\left(r,\mu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(r\right)\right)\mathrm{d}r-\mathfrak{t}{\int }_0^1{\displaystyle \frac{{(1-r)}^{{\ell }_1^{*}+{\ell }_2^{*}-1}}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*})}}M\left(r,\mu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(r\right)\right)\mathrm{d}r}\hfill \\ {\displaystyle +{\displaystyle \frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2^{*}+1}}{\mathsf{\Gamma }({\ell }_2^{*}+2)}}{\int }_0^1\frac{{(1-r)}^{{\ell }_1^{*}-1}}{\mathsf{\Gamma }\left({\ell }_1^{*}\right)}M\left(r,\mu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(r\right)\right)\mathrm{d}r,(\forall \mathfrak{t}\in \mathcal{O}).}\hfill \end{array}\right.$$

(10)

Based on Proposition 3 and because of our further proofs, we here define $H_1:\mathfrak{B}\to \mathfrak{B}$ and $H_2:\mathfrak{B}\to \mathfrak{B}$ by the following:

$$\begin{array}{cc}\hfill \left(H_1\nu \right)\left(\mathfrak{t}\right)& {\displaystyle ={\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)\mathrm{d}r}\hfill \\ \hfill & {\displaystyle -\mathfrak{t}{\int }_0^1{\displaystyle \frac{{(1-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)\mathrm{d}r}\hfill \\ \hfill & {\displaystyle +{\displaystyle \frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+2)}}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)\mathrm{d}r,}\hfill \end{array}$$

(11)

and

$$\begin{array}{cc}\hfill \left(H_2\mu \right)\left(\mathfrak{t}\right)& {\displaystyle ={\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1^{*}+{\ell }_2^{*}-1}}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*})}}M\left(r,\mu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(r\right)\right)\mathrm{d}r}\hfill \\ \hfill & {\displaystyle -\mathfrak{t}{\int }_0^1{\displaystyle \frac{{(1-r)}^{{\ell }_1^{*}+{\ell }_2^{*}-1}}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*})}}M\left(r,\mu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(r\right)\right)\mathrm{d}r}\hfill \\ \hfill & {\displaystyle +{\displaystyle \frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2^{*}+1}}{\mathsf{\Gamma }({\ell }_2^{*}+2)}}{\int }_0^1\frac{{(1-r)}^{{\ell }_1^{*}-1}}{\mathsf{\Gamma }\left({\ell }_1^{*}\right)}M\left(r,\mu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(r\right)\right)\mathrm{d}r.}\hfill \end{array}$$

(12)

Thus, the following system of operator equations is produced as follows.

$$\left\{\begin{array}{c}\mu \left(\mathfrak{t}\right)=H_1\nu \left(\mathfrak{t}\right),\hfill \\ \nu \left(\mathfrak{t}\right)=H_2\mu \left(\mathfrak{t}\right).\hfill \end{array}\right.$$

(13)

Now, we introduce $H^{🟉}:\mathfrak{B}\times \mathfrak{B}⟶\mathfrak{B}\times \mathfrak{B}$ defined by $H^{🟉}(\mu ,\nu )=(H_1\nu ,H_2\mu )$, and we find its fixed points. To perform this, we try to use the theory of fixed points due to the Banach and Krasnoselskii. Let us now remember these two fixed point theorems.

Theorem 1

(Banach, [35]). Let $\mathfrak{B}$ be a complete metric space and $\mathbb{H}:\mathfrak{B}\to \mathfrak{B}$ be a contraction. Then, there is a unique $\mu \in \mathfrak{B}$ such that $\mathbb{H}\left(\mu \right)=\mu .$

Theorem 2

(Krasnoselskii, [35]). Assume that $\mathbb{E}\ne \varnothing$ is a bounded, closed, and convex set in the Banach space $\mathfrak{B}$ and ${\mathcal{H}}_1,{\mathcal{H}}_2:\mathbb{E}\to \mathfrak{B}$ are two operators satisfying the following: (a) ${\mathcal{H}}_1{\mu }_1+{\mathcal{H}}_2{\mu }_2\in \mathbb{E},$ for each ${\mu }_1,{\mu }_2\in \mathbb{E}$; (b) ${\mathcal{H}}_1$ is a contraction; (c) ${\mathcal{H}}_2$ is compact and continuous. Then, there is $\mu \in \mathbb{E}$ such that ${\mathcal{H}}_1\mu +{\mathcal{H}}_2\mu =\mu .$

To continue the proof of our main results, we focus on the following hypotheses that are required.

Hypothese 1 (H1).

$M,K\in C(\mathcal{O}\times {\mathbb{R}}^{\ge 0}\times {\mathbb{R}}^{\ge 0},{\mathbb{R}}^{\ge 0})$for$\mathfrak{t}\in \mathcal{O}$and$({\mu }_1,{\mu }_2),({\widehat{\mu }}_1,{\widehat{\mu }}_2),({\nu }_1,{\nu }_2),$$({\widehat{\nu }}_1,{\widehat{\nu }}_2)\in {\mathbb{R}}^2.$

Hypothese 2 (H2).

There is$0<{\mathbf{C}}_K\in \mathbb{R}$satisfying the following.

$$|K(\mathfrak{t},{\mu }_1,{\mu }_2)-K(\mathfrak{t},\widehat{{\mu }_1},\widehat{{\mu }_2})|\le {\mathbf{C}}_K\left(|{\mu }_1-\widehat{{\mu }_1}|+|{\mu }_2-\widehat{{\mu }_2}|\right),\mathfrak{t}\in \mathcal{O},({\mu }_1,{\mu }_2),(\widehat{{\mu }_1},\widehat{{\mu }_2})\in {\mathbb{R}}^2.$$

Hypothese 3 (H3).

There is$0<{\mathbf{C}}_M\in \mathbb{R}$satisfying the following.

$$|M(\mathfrak{t},{\nu }_1,{\nu }_2)-K(\mathfrak{t},\widehat{{\nu }_1},\widehat{{\nu }_2})|\le {\mathbf{C}}_M\left(|{\nu }_1-\widehat{{\nu }_1}|+|{\nu }_2-\widehat{{\nu }_2}|\right),$$

for$\mathfrak{t}\in \mathcal{O}$and$({\nu }_1,{\nu }_2),(\widehat{{\nu }_1},\widehat{{\nu }_2})\in {\mathbb{R}}^2.$

Now, for more convenience in computations, we set the following.

$${\mathsf{\Delta }}_1=\frac{2}{\mathsf{\Gamma }({\ell }_1+{\ell }_2+1)}+\frac{1}{\mathsf{\Gamma }({\ell }_1+1)\mathsf{\Gamma }({\ell }_2+2)}+\frac{1}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }({\ell }_1+{\ell }_2-1)}+\frac{1}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }\left({\ell }_2\right)\mathsf{\Gamma }({\ell }_1+1)},$$

$${\mathsf{\Delta }}_2=\frac{2}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}+1)}+\frac{1}{\mathsf{\Gamma }({\ell }_1^{*}+1)\mathsf{\Gamma }({\ell }_2^{*}+2)}+\frac{1}{\mathsf{\Gamma }(3-{\ell }_2^{*})\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}-1)}+\frac{1}{\mathsf{\Gamma }(3-{\ell }_2^{*})\mathsf{\Gamma }\left({\ell }_2^{*}\right)\mathsf{\Gamma }({\ell }_1^{*}+1)}.$$

At this time, we can establish the following existence theorem.

Theorem 3.

Suppose that Hypotheses ($H1$)–($H3$) are valid together with ${\mathbf{C}}_K{\mathsf{\Delta }}_1<1$ and ${\mathbf{C}}_M{\mathsf{\Delta }}_2<1$. Then, the coupled sequential FBVPs of the Navier model (1) possess a unique solution.

Proof. 

First, we choose $\varrho \ge max\left\{\frac{{\mathsf{\Delta }}_1{\theta }_1}{1-{\mathbf{C}}_K{\mathsf{\Delta }}_1},\frac{{\mathsf{\Delta }}_2{\theta }_2}{1-{\mathbf{C}}_M{\mathsf{\Delta }}_2}\right\},$ where ${\mathsf{\Theta }}_1={max}_{\mathfrak{t}\in \mathcal{O}}|K(\mathfrak{t},0,0)|$ and ${\mathsf{\Theta }}_2={max}_{\mathfrak{t}\in \mathcal{O}}|M(\mathfrak{t},0,0)|$. We consider the bounded, closed, and convex set $\mathsf{\Lambda }$ defined by $\mathsf{\Lambda }=\left\{(\mu ,\nu )\in \mathfrak{B}\times \mathfrak{B}:\parallel (\mu ,\nu ){\parallel }_{\mathfrak{B}\times \mathfrak{B}}\le \varrho \right\}.$ From Hypotheses (H1) and (H2), we can write, for each $\mathfrak{t}\in \mathcal{O}$ and $(\mu ,\nu )\in \mathsf{\Lambda },$ that the following is the case.

$$\begin{array}{ccc}\hfill |H_1\nu \left(\mathfrak{t}\right)|& \le & \frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}|K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))|\mathrm{d}r\hfill \\ & & +\frac{\mathfrak{t}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}{\int }_0^1{(1-r)}^{{\ell }_1+{\ell }_2-1}|K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))|\mathrm{d}r\hfill \\ & & +\frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }\left({\ell }_1\right)\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1{(1-r)}^{{\ell }_1-1}|K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))|\mathrm{d}r\hfill \\ & \le & \frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}\left(|K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))-K(r,0,0)|+|K(r,0,0)|\right)\mathrm{d}r\hfill \\ & & +\frac{\mathfrak{t}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}{\int }_0^1{(1-r)}^{{\ell }_1+{\ell }_2-1}\left(|K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))-K(r,0,0)|+|K(r,0,0)|\right)\mathrm{d}r\hfill \\ & & +\frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }\left({\ell }_1\right)\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1{(1-r)}^{{\ell }_1-1}\left(|K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))-K(r,0,0)|+|K(r,0,0)|\right)\mathrm{d}r\hfill \\ & \le & \frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}\left({\mathbf{C}}_K{\parallel \nu \parallel }_{\mathfrak{B}}+{\mathsf{\Theta }}_1\right)\mathrm{d}r\hfill \\ & & +\frac{\mathfrak{t}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}{\int }_0^1{(1-r)}^{{\ell }_1+{\ell }_2-1}\left({\mathbf{C}}_K{\parallel \nu \parallel }_{\mathfrak{B}}+{\mathsf{\Theta }}_1\right)\mathrm{d}r\hfill \\ & & +\frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }\left({\ell }_1\right)\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1{(1-r)}^{{\ell }_1-1}\left({\mathbf{C}}_K{\parallel \nu \parallel }_{\mathfrak{B}}+{\mathsf{\Theta }}_1\right)\mathrm{d}r\hfill \\ & \le & \left({\mathbf{C}}_K\varrho +{\mathsf{\Theta }}_1\right)\left[\frac{2}{\mathsf{\Gamma }({\ell }_1+{\ell }_2+1)}+\frac{1}{\mathsf{\Gamma }({\ell }_1+1)\mathsf{\Gamma }({\ell }_2+2)}\right].\hfill \end{array}$$

(14)

By some simple calculations, we obtain the following.

$$\begin{array}{ccc}\hfill \left|{\left(H_1\nu \right)}^{″}\left(\mathfrak{t}\right)\right|& =& |\frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2-2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-3}K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))\mathrm{d}r\hfill \\ & & -\frac{{\ell }_2{\mathfrak{t}}^{{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_2+1)\mathsf{\Gamma }\left({\ell }_1\right)}{\int }_0^1{(1-r)}^{{\ell }_1-1}K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))\mathrm{d}r|\hfill \\ & \le & \frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2-2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-3}|K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))|\mathrm{d}r\hfill \\ & & +\frac{{\ell }_2{\mathfrak{t}}^{{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_2+1)\mathsf{\Gamma }\left({\ell }_1\right)}{\int }_0^1{(1-r)}^{{\ell }_1-1}|K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))|\mathrm{d}r\hfill \\ & \le & \left({\mathbf{C}}_K\varrho +{\mathsf{\Theta }}_1\right)\left[\frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2-1)}+\frac{1}{\mathsf{\Gamma }\left({\ell }_2\right)\mathsf{\Gamma }({\ell }_1+1)}\right].\hfill \end{array}$$

Then, by Definition 3, we obtain the following.

$$\begin{array}{ccc}\hfill \left|{}^C{\mathfrak{D}}^{{\ell }_2}\left(H_1\nu \right)\left(\mathfrak{t}\right)\right|& =& \left|\frac{1}{\mathsf{\Gamma }(2-{\ell }_2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{1-{\ell }_2}{\left(H_1\nu \right)}^{″}\left(r\right)\mathrm{d}r\right|\hfill \\ & \le & \frac{1}{\mathsf{\Gamma }(2-{\ell }_2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{1-{\ell }_2}|{\left(H_1\nu \right)}^{″}\left(r\right)|\mathrm{d}r\hfill \\ & \le & \left({\mathbf{C}}_K\varrho +{\mathsf{\Theta }}_1\right){\mathfrak{t}}^{2-{\ell }_2}\left[\frac{1}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }({\ell }_1+{\ell }_2-1)}+\frac{1}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }\left({\ell }_2\right)\mathsf{\Gamma }({\ell }_1+1)}\right]\hfill \\ & \le & \left({\mathbf{C}}_K\varrho +{\mathsf{\Theta }}_1\right)\left[\frac{1}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }({\ell }_1+{\ell }_2-1)}+\frac{1}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }\left({\ell }_2\right)\mathsf{\Gamma }({\ell }_1+1)}\right].\hfill \end{array}$$

(15)

Thus, Equations (14) and (15) provide the following.

$$\begin{array}{ccc}\hfill \parallel H_1{\nu \parallel }_{\mathfrak{B}}& \le & \left({\mathbf{C}}_K\varrho +{\mathsf{\Theta }}_1\right)[\frac{2}{\mathsf{\Gamma }({\ell }_1+{\ell }_2+1)}+\frac{1}{\mathsf{\Gamma }({\ell }_1+1)\mathsf{\Gamma }({\ell }_2+2)}\hfill \\ & & +\frac{1}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }({\ell }_1+{\ell }_2-1)}+\frac{1}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }\left({\ell }_2\right)\mathsf{\Gamma }({\ell }_1+1)}]\hfill \\ & \le & \left({\mathbf{C}}_K\varrho +{\mathsf{\Theta }}_1\right){\mathsf{\Delta }}_1\le \varrho .\hfill \end{array}$$

(16)

By following the same procedure, we can obtain the following.

$$\begin{array}{ccc}\hfill \parallel H_2{\mu \parallel }_{\mathfrak{B}}& \le & \left({\mathbf{C}}_M\varrho +{\mathsf{\Theta }}_2\right)[\frac{2}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}+1)}+\frac{1}{\mathsf{\Gamma }({\ell }_1^{*}+1)\mathsf{\Gamma }({\ell }_2^{*}+2)}\hfill \\ & & +\frac{1}{\mathsf{\Gamma }(3-{\ell }_2^{*})\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}-1)}+\frac{1}{\mathsf{\Gamma }(3-{\ell }_2^{*})\mathsf{\Gamma }\left({\ell }_2^{*}\right)\mathsf{\Gamma }({\ell }_1^{*}+1)}]\hfill \\ & \le & \left({\mathbf{C}}_M\varrho +{\mathsf{\Theta }}_2\right){\mathsf{\Delta }}_2\le \varrho .\hfill \end{array}$$

(17)

Consequently, Equations (16) and (17) result in $\parallel H^{🟉}(\mu ,\nu ){\parallel }_{\mathfrak{B}\times \mathfrak{B}}\le \varrho .$ Therefore $H^{🟉}\left(\mathsf{\Lambda }\right)\subseteq \mathsf{\Lambda }$.

We prove that $H^{🟉}$ is a contraction. For this aim, let us consider $(\mu ,\nu ),(\widehat{\mu },\widehat{\nu })\in \mathsf{\Lambda }$, and $\mathfrak{t}\in \mathcal{O}$ arbitrarily. Then, we have the following.

$$\begin{array}{ccc}\hfill & & |H_1\nu \left(\mathfrak{t}\right)-H_1\widehat{\nu }\left(\mathfrak{t}\right)|\hfill \\ & \le & \frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}|K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))-K(r,\widehat{\nu }\left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\widehat{\nu }\left(r\right))|\mathrm{d}r\hfill \\ & & +\frac{\mathfrak{t}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}{\int }_0^1{(1-r)}^{{\ell }_1+{\ell }_2-1}|K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))-K(r,\widehat{\nu }\left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\widehat{\nu }\left(r\right))|\mathrm{d}r\hfill \\ & & +\frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }\left({\ell }_1\right)\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1{(1-r)}^{{\ell }_1-1}|K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))-K(r,\widehat{\nu }\left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\widehat{\nu }\left(r\right))|\mathrm{d}r\hfill \\ & \le & \frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}{\mathbf{C}}_K\left(|\nu \left(r\right)-\widehat{\nu }\left(r\right)|+|{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)-{}^C{\mathfrak{D}}^{{\ell }_2}\widehat{\nu }\left(r\right)|\right)\mathrm{d}r\hfill \\ & & +\frac{\mathfrak{t}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}{\int }_0^1{(1-r)}^{{\ell }_1+{\ell }_2-1}{\mathbf{C}}_K\left(|\nu \left(r\right)-\widehat{\nu }\left(r\right)|+|{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)-{}^C{\mathfrak{D}}^{{\ell }_2}\widehat{\nu }\left(r\right)|\right)\mathrm{d}r\hfill \\ & & +\frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }\left({\ell }_1\right)\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1{(1-r)}^{{\ell }_1-1}{\mathbf{C}}_K\left(|\nu \left(r\right)-\widehat{\nu }\left(r\right)|+|{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)-{}^C{\mathfrak{D}}^{{\ell }_2}\widehat{\nu }\left(r\right)|\right)\mathrm{d}r\hfill \\ & \le & {\mathbf{C}}_K\left[\frac{2}{\mathsf{\Gamma }({\ell }_1+{\ell }_2+1)}+\frac{1}{\mathsf{\Gamma }({\ell }_2+2)\mathsf{\Gamma }({\ell }_1+1)}\right]\parallel \nu -\widehat{\nu }{\parallel }_{\mathfrak{B}}.\hfill \end{array}$$

(18)

On the other hand, since the following is the case:

$$\begin{array}{ccc}\hfill & & |{\left(H_1\nu \right)}^{″}\left(\mathfrak{t}\right)-{\left(H_1\widehat{\nu }\right)}^{″}\left(\mathfrak{t}\right)|\hfill \\ & \le & \frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2-2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-3}|K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))-K(r,\widehat{\nu }\left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\widehat{\nu }\left(r\right))|\mathrm{d}r\hfill \\ & & +\frac{{\ell }_2{\mathfrak{t}}^{{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_2+1)\mathsf{\Gamma }\left({\ell }_1\right)}{\int }_0^1{(1-r)}^{{\ell }_1-1}|K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))-K(r,\widehat{\nu }\left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\widehat{\nu }\left(r\right))|\mathrm{d}r\hfill \\ & \le & \frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2-2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-3}{\mathbf{C}}_K\left(|\nu \left(r\right)-\widehat{\nu }\left(r\right)|+|{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)-{}^C{\mathfrak{D}}^{{\ell }_2}\widehat{\nu }\left(r\right)|\right)\mathrm{d}r\hfill \\ & & +\frac{{\ell }_2{\mathfrak{t}}^{{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_2+1)\mathsf{\Gamma }\left({\ell }_1\right)}{\int }_0^1{(1-r)}^{{\ell }_1-1}{\mathbf{C}}_K\left(|\nu \left(r\right)-\widehat{\nu }\left(r\right)|+|{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)-{}^C{\mathfrak{D}}^{{\ell }_2}\widehat{\nu }\left(r\right)|\right)\mathrm{d}r\hfill \\ & \le & {\mathbf{C}}_K\left[\frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2-1)}+\frac{1}{\mathsf{\Gamma }\left({\ell }_2\right)\mathsf{\Gamma }({\ell }_1+1)}\right]\parallel \nu -\widehat{\nu }{\parallel }_{\mathfrak{B}},\hfill \end{array}$$

(19)

by using (19), we obtain the following.

$$\begin{array}{ccc}\hfill & & |{}^C{\mathfrak{D}}^{{\ell }_2}\left(H_1\nu \right)\left(\mathfrak{t}\right)-{}^C{\mathfrak{D}}^{{\ell }_2}\left(H_1\widehat{\nu }\right)\left(\mathfrak{t}\right)|\hfill \\ & =& |\frac{1}{\mathsf{\Gamma }(2-{\ell }_2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{1-{\ell }_2}{\left(H_1\nu \right)}^{″}\left(r\right)\mathrm{d}r-\frac{1}{\mathsf{\Gamma }(2-{\ell }_2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{1-{\ell }_2}{\left(H_1\widehat{\nu }\right)}^{″}\left(r\right)\mathrm{d}r|\hfill \\ & \le & \frac{1}{\mathsf{\Gamma }(2-{\ell }_2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{1-{\ell }_2}|{\left(H_1\nu \right)}^{″}\left(r\right)-{\left(H_1\widehat{\nu }\right)}^{″}\left(r\right)|\mathrm{d}r\hfill \\ & \le & {\mathbf{C}}_K\left[\frac{1}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }({\ell }_1+{\ell }_2-1)}+\frac{1}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }\left({\ell }_2\right)\mathsf{\Gamma }({\ell }_1+1)}\right]\parallel \nu -\widehat{\nu }{\parallel }_{\mathfrak{B}}.\hfill \end{array}$$

(20)

From Equations (18) and (20), we obtain the following.

$$\parallel H_1\nu -H_1\widehat{\nu }{\parallel }_{\mathfrak{B}}\le {\mathbf{C}}_K{\mathsf{\Delta }}_1\parallel \nu -\widehat{\nu }{\parallel }_{\mathfrak{B}}.$$

(21)

By the same arguments, we obtain the following.

$$\parallel H_2\mu -H_2\widehat{\mu }{\parallel }_{\mathfrak{B}}\le {\mathbf{C}}_M{\mathsf{\Delta }}_2\parallel \mu -\widehat{\mu }{\parallel }_{\mathfrak{B}}.$$

(22)

Therefore, Equations (21) and (22) together with the assumptions ${\mathbf{C}}_K{\mathsf{\Delta }}_1<1$ and ${\mathbf{C}}_M{\mathsf{\Delta }}_2<1$ result in $\parallel H^{🟉}(\mu ,\nu )-H^{🟉}(\widehat{\mu },\widehat{\nu }){\parallel }_{\mathfrak{B}\times \mathfrak{B}}\le \parallel (\mu ,\nu )-(\widehat{\mu },\widehat{\nu }){\parallel }_{\mathfrak{B}\times \mathfrak{B}}.$ This means that $H^{🟉}$ is a contraction. Hence, Theorem 1 ensured that the coupled sequential FBVPs of Navier model (1) possess a unique solution; thus, the proof is now completed. □

In order to continue the establishment of the main existence result, we need to decompose the operators $H_1$ and $H_2$ as follows: $H_1={\widehat{\mathsf{\Phi }}}_1+{\widehat{\mathsf{\Psi }}}_1$ and $H_2={\widehat{\mathsf{\Phi }}}_2+{\widehat{\mathsf{\Psi }}}_2$, where the following is the case.

$$\left\{\begin{array}{c}{\displaystyle {\widehat{\mathsf{\Phi }}}_1\nu \left(\mathfrak{t}\right)={\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)\mathrm{d}r,}\hfill \\ {\displaystyle {\widehat{\mathsf{\Psi }}}_1\nu \left(\mathfrak{t}\right)=-\mathfrak{t}{\int }_0^1{\displaystyle \frac{{(1-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)\mathrm{d}r}\hfill \\ {\displaystyle +{\displaystyle \frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+2)}}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)\mathrm{d}r,}\hfill \\ {\displaystyle {\widehat{\mathsf{\Phi }}}_2\mu \left(\mathfrak{t}\right)={\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1^{*}+{\ell }_2^{*}-1}}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*})}}M\left(r,\mu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(r\right)\right)\mathrm{d}r,}\hfill \\ {\displaystyle {\widehat{\mathsf{\Psi }}}_2\mu \left(\mathfrak{t}\right)=-\mathfrak{t}{\int }_0^1{\displaystyle \frac{{(1-r)}^{{\ell }_1^{*}+{\ell }_2^{*}-1}}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*})}}M\left(r,\mu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(r\right)\right)\mathrm{d}r}\hfill \\ {\displaystyle +{\displaystyle \frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2^{*}+1}}{\mathsf{\Gamma }({\ell }_2^{*}+2)}}{\int }_0^1\frac{{(1-r)}^{{\ell }_1^{*}-1}}{\mathsf{\Gamma }\left({\ell }_1^{*}\right)}M\left(r,\mu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(r\right)\right)\mathrm{d}r.}\hfill \end{array}\right.$$

(23)

In this case, the operator $H^{🟉}$ can be expressed as $H^{🟉}=\widehat{\mathsf{\Phi }}+\widehat{\mathsf{\Psi }},$ where $\widehat{\mathsf{\Phi }}$ and $\widehat{\mathsf{\Psi }}$ are given by the following.

$$\widehat{\mathsf{\Phi }}(\mu ,\nu )=({\widehat{\mathsf{\Phi }}}_1\nu ,{\widehat{\mathsf{\Phi }}}_2\mu )\mathrm{and}\widehat{\mathsf{\Psi }}(\mu ,\nu )=({\widehat{\mathsf{\Psi }}}_1\nu ,{\widehat{\mathsf{\Psi }}}_2\mu ).$$

In addition, suppose that the functions K and M satisfy the following assumption.

Hypothese 4 (H4).

There exists${{\rm Y}}_K,{{\rm Y}}_M,{\Omega }_K,{\Omega }_M\in {\mathbb{R}}^{+}$such that the following is the case:

$$|K(\mathfrak{t},\nu \left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(\mathfrak{t}\right))|\le {{\rm Y}}_K\left(|\nu \left(\mathfrak{t}\right)|+|{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(\mathfrak{t}\right)|\right)+{\Omega }_K,$$

and

$$|M(\mathfrak{t},\mu \left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(\mathfrak{t}\right))|\le {{\rm Y}}_M\left(|\mu \left(\mathfrak{t}\right)|+|{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(\mathfrak{t}\right)|\right)+{\Omega }_M,$$

for all$\mathfrak{t}\in \mathcal{O}$and for each$(\mu ,\nu )\in {\mathfrak{B}}^2.$

Theorem 4.

Suppose that the assumptions Hypotheses ($H1$)–($H4$) hold. Furthermore, the conditions

$$\begin{array}{c}\hfill {\mathbf{C}}_K\left[\frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2+1)}+\frac{1}{\mathsf{\Gamma }({\ell }_2+2)\mathsf{\Gamma }({\ell }_1+1)}+\frac{1}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }\left({\ell }_2\right)\mathsf{\Gamma }({\ell }_1+1)}\right]<1,\end{array}$$

(24)

and

$$\begin{array}{c}\hfill {\mathbf{C}}_M\left[\frac{1}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}+1)}+\frac{1}{\mathsf{\Gamma }({\ell }_2^{*}+2)\mathsf{\Gamma }({\ell }_1^{*}+1)}+\frac{1}{\mathsf{\Gamma }(3-{\ell }_2^{*})\mathsf{\Gamma }\left({\ell }_2^{*}\right)\mathsf{\Gamma }({\ell }_1^{*}+1)}\right]<1,\end{array}$$

(25)

are valid. Then, the coupled sequential FBVPs of Navier model (1) has at least one solution.

Proof. 

First of all, the continuity property of the operator $H^{🟉}$ is obtained from that of the functions K and M. Let the set D be bounded in $\mathsf{\Lambda }\in \mathfrak{B}\times \mathfrak{B}$. Then, in view of the assumption (H4), we have the following for any $\mathfrak{t}\in \mathcal{O}$ and all $(\mu ,\nu )\in D$.

$$\begin{array}{cc}\hfill |{\widehat{\mathsf{\Phi }}}_1\nu \left(\mathfrak{t}\right)|\le & {\displaystyle {\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}|K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)|\mathrm{d}r}\hfill \\ \hfill \le & {\displaystyle \frac{{{\rm Y}}_K\parallel \nu {\parallel }_{\mathfrak{B}}+{\Omega }_K}{\mathsf{\Gamma }({\ell }_1+{\ell }_2+1)}.}\hfill \end{array}$$

(26)

On the other hand, by

$$\begin{array}{cc}\hfill \left|{\left({\widehat{\mathsf{\Phi }}}_1\nu \right)}^{″}\left(\mathfrak{t}\right)\right|\le & {\displaystyle {\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-3}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2-2)}}|K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)|\mathrm{d}r}\hfill \\ \hfill \le & {\displaystyle \frac{{{\rm Y}}_K\parallel \nu {\parallel }_{\mathfrak{B}}+{\Omega }_K}{\mathsf{\Gamma }({\ell }_1+{\ell }_2-1)},}\hfill \end{array}$$

and by Definition 3, we obtain the following.

$$\begin{array}{cc}\hfill \left|{}^C{\mathfrak{D}}^{{\ell }_2}\left({\widehat{\mathsf{\Phi }}}_1\nu \right)\left(\mathfrak{t}\right)\right|\le & {\displaystyle {\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{1-{\ell }_2}}{\mathsf{\Gamma }(2-{\ell }_2)}}\left|{\left({\widehat{\mathsf{\Phi }}}_1\nu \right)}^{″}\left(r\right)\right|\mathrm{d}r}\hfill \\ \hfill \le & {\displaystyle \frac{{{\rm Y}}_K\parallel \nu {\parallel }_{\mathfrak{B}}+{\Omega }_K}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }({\ell }_1+{\ell }_2-1)}.}\hfill \end{array}$$

(27)

Thus, Equations (26) and (27) provide the following.

$$\begin{array}{c}\hfill \parallel {\widehat{\mathsf{\Phi }}}_1\nu {\parallel }_{\mathfrak{B}} \le \left({{\rm Y}}_K\parallel \nu {\parallel }_{\mathfrak{B}}+{\Omega }_K\right)\left[\frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2+1)}+\frac{1}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }({\ell }_1+{\ell }_2-1)}\right].\end{array}$$

(28)

In a similar manner, we have the following.

$$\begin{array}{c}\hfill \parallel {\widehat{\mathsf{\Phi }}}_2\mu {\parallel }_{\mathfrak{B}} \le \left({{\rm Y}}_M\parallel \mu {\parallel }_{\mathfrak{B}}+{\Omega }_M\right)\left[\frac{1}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}+1)}+\frac{1}{\mathsf{\Gamma }(3-{\ell }_2^{*})\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}-1)}\right].\end{array}$$

(29)

Therefore, Equations (28) and (29) provide the boundedness of $\widehat{\mathsf{\Phi }}\left(D\right).$

In the sequel, we shall investigate that $\widehat{\mathsf{\Phi }}$ is equicontinuous. Let ${\tau }_1,{\tau }_2\in \mathcal{O}$ with ${\tau }_1<{\tau }_2$ and $(\mu ,\nu )\in \mathfrak{B}\times \mathfrak{B}.$ We have the following.

$$\begin{array}{cc}\hfill |{\widehat{\mathsf{\Phi }}}_1\nu \left({\tau }_1\right)-{\widehat{\mathsf{\Phi }}}_1\nu \left({\tau }_2\right)|\le & {\displaystyle {\displaystyle \frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}[{\int }_0^{{\tau }_1}\left[{({\tau }_2-r)}^{{\ell }_1+{\ell }_2-1}-{({\tau }_1-r)}^{{\ell }_1+{\ell }_2-1}\right]|K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)|\mathrm{d}r}\hfill \\ \hfill & +{\int }_{{\tau }_1}^{{\tau }_2}{({\tau }_2-r)}^{{\ell }_1+{\ell }_2-1}|K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)|\mathrm{d}r]\hfill \\ \hfill \le & {\displaystyle \frac{{{\rm Y}}_K\parallel \nu {\parallel }_{\mathfrak{B}}+{\Omega }_K}{\mathsf{\Gamma }({\ell }_1+{\ell }_2+1)}\left({\tau }_2^{{\ell }_1+{\ell }_2}-{\tau }_1^{{\ell }_1+{\ell }_2}\right).}\hfill \end{array}$$

(30)

In a similar manner, we obtain the following.

$$\begin{array}{c}\hfill {\displaystyle |{\widehat{\mathsf{\Phi }}}_2\mu \left({\tau }_1\right)-{\widehat{\mathsf{\Phi }}}_2\mu \left({\tau }_2\right)|\le \frac{{{\rm Y}}_M\parallel \mu {\parallel }_{\mathfrak{B}}+{\Omega }_M}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}+1)}\left({\tau }_2^{{\ell }_1^{*}+{\ell }_2^{*}}-{\tau }_1^{{\ell }_1^{*}+{\ell }_2^{*}}\right).}\end{array}$$

(31)

We also have the following.

$$\begin{array}{cc}\hfill & |{}^C{\mathfrak{D}}^{{\ell }_2}\left({\widehat{\mathsf{\Phi }}}_1\nu \right)\left({\tau }_1\right)-{}^C{\mathfrak{D}}^{{\ell }_2}\left({\widehat{\mathsf{\Phi }}}_1\nu \right)\left({\tau }_2\right)|\hfill \\ \hfill \le & {\displaystyle {\displaystyle \frac{1}{\mathsf{\Gamma }(2-{\ell }_2)}}\left[{\int }_0^{{\tau }_1}\left[{({\tau }_1-r)}^{1-{\ell }_2}-{({\tau }_2-r)}^{1-{\ell }_2}\right]\left|{\left({\widehat{\mathsf{\Phi }}}_1\nu \right)}^{″}\left(r\right)\right|\mathrm{d}r+{\int }_{{\tau }_1}^{{\tau }_2}{({\tau }_2-r)}^{1-{\ell }_2}\left|{\left({\widehat{\mathsf{\Phi }}}_1\nu \right)}^{″}\left(r\right)\right|\mathrm{d}r\right]}\hfill \\ \hfill \le & {\displaystyle \frac{2\left({{\rm Y}}_K\parallel \nu {\parallel }_{\mathfrak{B}}+{\Omega }_K\right)}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }({\ell }_1+{\ell }_2-1)}{\left({\tau }_2-{\tau }_1\right)}^{2-{\ell }_2}.}\hfill \end{array}$$

(32)

By using the same arguments, we obtain the following.

$${\displaystyle |{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\left({\widehat{\mathsf{\Phi }}}_2\mu \right)\left({\tau }_1\right)-{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\left({\widehat{\mathsf{\Phi }}}_2\mu \right)\left({\tau }_2\right)|\le \frac{2\left({{\rm Y}}_M\parallel \mu {\parallel }_{\mathfrak{B}}+{\Omega }_M\right)}{\mathsf{\Gamma }(3-{\ell }_2^{*})\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}-1)}{\left({\tau }_2-{\tau }_1\right)}^{2-{\ell }_2^{*}}.}$$

(33)

From inequalities (Equations (30)–(33)), we conclude

$$\begin{array}{c}\hfill \parallel {\widehat{\mathsf{\Phi }}}_1\nu \left({\tau }_1\right)-{\widehat{\mathsf{\Phi }}}_1\nu \left({\tau }_2\right){\parallel }_{\mathfrak{B}} \le \left({{\rm Y}}_K\parallel \nu {\parallel }_{\mathfrak{B}}+{\Omega }_K\right)\left[\frac{{\tau }_2^{{\ell }_1+{\ell }_2}-{\tau }_1^{{\ell }_1+{\ell }_2}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2+1)}+\frac{2{({\tau }_2-{\tau }_1)}^{2-{\ell }_2}}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }({\ell }_1+{\ell }_2-1)}\right],\end{array}$$

(34)

and

$$\begin{array}{c}\hfill \parallel {\widehat{\mathsf{\Phi }}}_2\mu \left({\tau }_1\right)-{\widehat{\mathsf{\Phi }}}_2\mu \left({\tau }_2\right){\parallel }_{\mathfrak{B}} \le \left({{\rm Y}}_M\parallel \mu {\parallel }_{\mathfrak{B}}+{\Omega }_M\right)\left[\frac{{\tau }_2^{{\ell }_1^{*}+{\ell }_2^{*}}-{\tau }_1^{{\ell }_1^{*}+{\ell }_2^{*}}}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}+1)}+\frac{2{({\tau }_2-{\tau }_1)}^{2-{\ell }_2^{*}}}{\mathsf{\Gamma }(3-{\ell }_2^{*})\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}-1)}\right].\end{array}$$

(35)

The obtained inequalities (34)–(35) provide $\parallel {\widehat{\mathsf{\Phi }}}_1\nu \left({\tau }_1\right)-{\widehat{\mathsf{\Phi }}}_1\nu \left({\tau }_2\right){\parallel }_{\mathfrak{B}}\to 0$ and $\parallel {\widehat{\mathsf{\Phi }}}_2\mu \left({\tau }_1\right)-{\widehat{\mathsf{\Phi }}}_2\mu \left({\tau }_2\right){\parallel }_{\mathfrak{B}}\to 0$ if ${\tau }_1-{\tau }_2\to 0.$ Therefore, by Arzela–Ascoli’s theorem, it follows that $\widehat{\mathsf{\Phi }}$ is continuous and compact.

Finally, we prove the last part of the theorem regarding the operator $\widehat{\mathsf{\Psi }}$. For each $(\nu ,\widehat{\nu })\in \mathfrak{B}\times \mathfrak{B},$ we have the following.

$$\begin{array}{cc}\hfill |{\widehat{\mathsf{\Psi }}}_1\nu \left(\mathfrak{t}\right)-{\widehat{\mathsf{\Psi }}}_1\widehat{\nu }\left(\mathfrak{t}\right)|& {\displaystyle \le {\displaystyle \frac{\mathfrak{t}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}{\int }_0^1{(1-r)}^{{\ell }_1+{\ell }_2-1}|K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)-K\left(r,\widehat{\nu }\left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\widehat{\nu }\left(r\right)\right)|\mathrm{d}r}\hfill \\ \hfill & +{\displaystyle \frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+1)\mathsf{\Gamma }\left({\ell }_1\right)}}{\int }_0^1{(1-r)}^{{\ell }_1-1}|K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)-K\left(r,\widehat{\nu }\left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\widehat{\nu }\left(r\right)\right)|\mathrm{d}r\hfill \\ \hfill & \le {\mathbf{C}}_K\left[{\displaystyle \frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2+1)}}+{\displaystyle \frac{1}{\mathsf{\Gamma }({\ell }_2+1)\mathsf{\Gamma }({\ell }_1+1)}}\right]\parallel \nu -\widehat{\nu }{\parallel }_{\mathfrak{B}}.\hfill \end{array}$$

(36)

Similarly, we have the following.

$$\begin{array}{c}\hfill |{\widehat{\mathsf{\Psi }}}_2\mu \left(\mathfrak{t}\right)-{\widehat{\mathsf{\Psi }}}_2\widehat{\mu }\left(\mathfrak{t}\right)|\le {\mathbf{C}}_M\left[{\displaystyle \frac{1}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}+1)}}+{\displaystyle \frac{1}{\mathsf{\Gamma }({\ell }_2^{*}+1)\mathsf{\Gamma }({\ell }_1^{*}+1)}}\right]\parallel \mu -\widehat{\mu }{\parallel }_{\mathfrak{B}}.\end{array}$$

(37)

We also have the following:

$$\begin{array}{cc}\hfill |{\left({\widehat{\mathsf{\Psi }}}_1\nu \right)}^{″}\left(\mathfrak{t}\right)-{\left({\widehat{\mathsf{\Psi }}}_1\widehat{\nu }\right)}^{″}\left(\mathfrak{t}\right)|\le & {\displaystyle \frac{{\mathbf{C}}_K}{\mathsf{\Gamma }\left({\ell }_2\right)\mathsf{\Gamma }({\ell }_1+1)}}\parallel \nu -\widehat{\nu }{\parallel }_{\mathfrak{B}},\hfill \end{array}$$

(38)

and

$$\begin{array}{c}\hfill |{\left({\widehat{\mathsf{\Psi }}}_2\mu \right)}^{″}\left(\mathfrak{t}\right)-{\left({\widehat{\mathsf{\Psi }}}_2\widehat{\mu }\right)}^{″}\left(\mathfrak{t}\right)|\le {\displaystyle \frac{{\mathbf{C}}_M}{\mathsf{\Gamma }\left({\ell }_2^{*}\right)\mathsf{\Gamma }({\ell }_1^{*}+1)}}\parallel \mu -\widehat{\mu }{\parallel }_{\mathfrak{B}}.\end{array}$$

(39)

Then, the following is the case:

$$\begin{array}{cc}\hfill |{}^C{\mathfrak{D}}^{{\ell }_2}{\widehat{\mathsf{\Psi }}}_1\nu \left(\mathfrak{t}\right)-{}^C{\mathfrak{D}}^{{\ell }_2}{\widehat{\mathsf{\Psi }}}_1\widehat{\nu }\left(\mathfrak{t}\right)|\le & {\displaystyle {\displaystyle \frac{1}{\mathsf{\Gamma }(2-{\ell }_2)}}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{1-{\ell }_2}|{\left({\widehat{\mathsf{\Psi }}}_1\nu \right)}^{″}\left(r\right)-{\left({\widehat{\mathsf{\Psi }}}_1\widehat{\nu }\right)}^{″}\left(r\right)|\mathrm{d}r}\hfill \\ \hfill \le & {\displaystyle \frac{{\mathbf{C}}_K}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }\left({\ell }_2\right)\mathsf{\Gamma }({\ell }_1+1)}}\parallel \nu -\widehat{\nu }{\parallel }_{\mathfrak{B}},\hfill \end{array}$$

(40)

and

$$\begin{array}{c}\hfill |{}^C{\mathfrak{D}}^{{\ell }_2^{*}}{\widehat{\mathsf{\Psi }}}_2\mu \left(\mathfrak{t}\right)-{}^C{\mathfrak{D}}^{{\ell }_2^{*}}{\widehat{\mathsf{\Psi }}}_2\widehat{\mu }\left(\mathfrak{t}\right)|\le {\displaystyle \frac{{\mathbf{C}}_M}{\mathsf{\Gamma }(3-{\ell }_2^{*})\mathsf{\Gamma }\left({\ell }_2^{*}\right)\mathsf{\Gamma }({\ell }_1^{*}+1)}}\parallel \mu -\widehat{\mu }{\parallel }_{\mathfrak{B}}.\end{array}$$

(41)

Consequently, in view of Equations (36), (37), (40), and (41), we obtain

$$|{\widehat{\mathsf{\Psi }}}_1\nu \left(\mathfrak{t}\right)-{\widehat{\mathsf{\Psi }}}_1\widehat{\nu }\left(\mathfrak{t}\right)|\le {\mathbf{C}}_K\left[{\displaystyle \frac{1}{\mathsf{\Gamma }({\ell }_1+{\ell }_2+1)}}+{\displaystyle \frac{1}{\mathsf{\Gamma }({\ell }_2+2)\mathsf{\Gamma }({\ell }_1+1)}}+{\displaystyle \frac{1}{\mathsf{\Gamma }(3-{\ell }_2)\mathsf{\Gamma }\left({\ell }_2\right)\mathsf{\Gamma }({\ell }_1+1)}}\right]\parallel \nu -\widehat{\nu }{\parallel }_{\mathfrak{B}},$$

(42)

and

$$|{\widehat{\mathsf{\Psi }}}_2\mu \left(\mathfrak{t}\right)-{\widehat{\mathsf{\Psi }}}_2\widehat{\mu }\left(\mathfrak{t}\right)|\le {\mathbf{C}}_M\left[{\displaystyle \frac{1}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}+1)}}+{\displaystyle \frac{1}{\mathsf{\Gamma }({\ell }_2^{*}+2)\mathsf{\Gamma }({\ell }_1^{*}+1)}}+{\displaystyle \frac{1}{\mathsf{\Gamma }(3-{\ell }_2^{*})\mathsf{\Gamma }\left({\ell }_2^{*}\right)\mathsf{\Gamma }({\ell }_1^{*}+1)}}\right]\parallel \mu -\widehat{\mu }{\parallel }_{\mathfrak{B}}.$$

(43)

Therefore, from Equations (42) and (43) together with the conditions (24) and (25), it follows that $\widehat{\mathsf{\Psi }}$ is a contraction. Hence, by the conclusion of Theorem 2, we conclude $H^{🟉}$ has at least one fixed point which is the solution of the coupled sequential FBVPs of Navier model (1). □

4. Results for H-U-Type Stability

In the present section, we investigate some sufficient conditions to obtain H-U-stability results of the solutions to the coupled sequential FBVPs of Navier model (1). Before starting work in this regard, we present some auxiliary definitions.

Definition 4

([36,37]). Let $H^{🟉}:\mathfrak{B}\to \mathfrak{B}$. Then, the operator equation of the following:

$$\begin{array}{c}\hfill H^{🟉}\mathfrak{p}=\mathfrak{p},\mathfrak{p}\in \mathfrak{B},\end{array}$$

(44)

is said to be H-U-stable, if for the given inequality

$$|\mathfrak{p}\left(\mathfrak{t}\right)-H^{🟉}\mathfrak{p}\left(\mathfrak{t}\right)|\le ϵ,\forall \mathfrak{t}\in \mathcal{O},$$

there is some $w_{H^{🟉}}>0$ such that for each solution $\mathfrak{p}\in C(\mathcal{O},\mathbb{R})$ of the Equation (44), we can find some $\widehat{h}\in C(\mathcal{O},\mathbb{R})$ satisfying (44) such that the following is the case.

$$|\mathfrak{p}\left(\mathfrak{t}\right)-\widehat{h}\left(\mathfrak{t}\right)|\le w_{H^{🟉}}ϵ,\forall \mathfrak{t}\in \mathcal{O}.$$

Now, based on above definition, we define two operators $H_i:\mathfrak{B}⟶\mathfrak{B},$$i\in \{1,2\}.$ The coupled system of the following:

$$\left\{\begin{array}{c}\mu \left(\mathfrak{t}\right)=H_1\nu \left(\mathfrak{t}\right),\hfill \\ \nu \left(\mathfrak{t}\right)=H_2\mu \left(\mathfrak{t}\right),\hfill \end{array}\right.$$

(45)

is called H-U-stable if, for the system, the following is the case.

$$\left\{\begin{array}{c}|\mu \left(\mathfrak{t}\right)-H_1\nu \left(\mathfrak{t}\right)|\le {ϵ}_1,\mathfrak{t}\in \mathcal{O},\hfill \\ |\nu \left(\mathfrak{t}\right)-H_2\mu \left(\mathfrak{t}\right)|\le {ϵ}_2,\mathfrak{t}\in \mathcal{O},\hfill \end{array}\right.$$

(46)

Two constants $w_{H_1},w_{H_2}>0$ can be found provided that for each $(\mu ,\nu )$ satisfying (45), a unique solution $(\widehat{h},\tilde{h})$ exists for the system (45) such that the following is the case.

$$\left\{\begin{array}{c}|\mu \left(\mathfrak{t}\right)-\widehat{h}\left(\mathfrak{t}\right)|\le w_{H_1}{ϵ}_1,\forall \mathfrak{t}\in \mathcal{O},\hfill \\ |\nu \left(\mathfrak{t}\right)-\tilde{h}\left(\mathfrak{t}\right)|\le w_{H_2}{ϵ}_2,\forall \mathfrak{t}\in \mathcal{O}.\hfill \end{array}\right.$$

Remark 1.

There are two functions $\phi ,\chi \in C(\mathcal{O},\mathbb{R})$ with respect to μ and ν, respectively, and that satisfy the following:

$$\left\{\begin{array}{c}|\phi \left(\mathfrak{t}\right)|\le {ϵ}_1,for\mathfrak{t}\in \mathcal{O},\hfill \\ |\chi \left(\mathfrak{t}\right)|\le {ϵ}_2,for\mathfrak{t}\in \mathcal{O},\hfill \end{array}\right.$$

and

$$\left\{\begin{array}{c}{}^C{\mathfrak{D}}^{{\ell }_1}\left({}^C{\mathfrak{D}}^{{\ell }_2}\mu \right)\left(\mathfrak{t}\right)=K\left(\mathfrak{t},\nu \left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(\mathfrak{t}\right)\right) + \phi \left(\mathfrak{t}\right),(\mathfrak{t}\in \mathcal{O}:=[0,1]),\hfill \\ {}^C{\mathfrak{D}}^{{\ell }_1^{*}}\left({}^C{\mathfrak{D}}^{{\ell }_2^{*}}\nu \right)\left(\mathfrak{t}\right)=M\left(\mathfrak{t},\mu \left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(\mathfrak{t}\right)\right) + \chi \left(\mathfrak{t}\right),(\mathfrak{t}\in \mathcal{O}:=[0,1]).\hfill \end{array}\right.$$

Lemma 1.

Suppose that $(\mu ,\nu )\in {\left(C(\mathcal{O},\mathbb{R})\right)}^2$ is a solution of (46) with the following.

$$\left\{\begin{array}{c}|{\mu }^{″}\left(\mathfrak{t}\right)-{\left(H_1\nu \right)}^{″}\left(\mathfrak{t}\right)|\le {\gamma }_1,{\gamma }_1\in {\mathbb{R}}^{>0},\mathfrak{t}\in \mathcal{O},\hfill \\ |{\nu }^{″}\left(\mathfrak{t}\right)-{\left(H_2\mu \right)}^{″}\left(\mathfrak{t}\right)|\le {\gamma }_2,{\gamma }_2\in {\mathbb{R}}^{>0},\mathfrak{t}\in \mathcal{O}.\hfill \end{array}\right.$$

Then, the following system of inequalities holds:

$$\left\{\begin{array}{c}\parallel \mu -H_1\nu {\parallel }_{\mathfrak{B}}\le {\eta }_1{ϵ}_1+{\displaystyle \frac{{\gamma }_1}{\mathsf{\Gamma }(3-{\ell }_2)}},\hfill \\ \parallel \nu -H_2\mu {\parallel }_{\mathfrak{B}}\le {\eta }_2{ϵ}_2+{\displaystyle \frac{{\gamma }_2}{\mathsf{\Gamma }(3-{\ell }_2^{*})}},\hfill \end{array}\right.$$

where ${\displaystyle {\eta }_1=\frac{2}{\mathsf{\Gamma }({\ell }_1+{\ell }_2+1)}+\frac{1}{\mathsf{\Gamma }({\ell }_2+2)\mathsf{\Gamma }({\ell }_1+1)},}$ ${\displaystyle {\eta }_2=\frac{2}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*}+1)}+\frac{1}{\mathsf{\Gamma }({\ell }_2^{*}+2)\mathsf{\Gamma }({\ell }_1^{*}+1)}}$, and $H_1$ and $H_2$ are defined by Equations (11) and (12), respectively.

Proof. 

In view of the condition $\left(2\right)$ (Remark 1), we have the following.

$$\left\{\begin{array}{c}{}^C{\mathfrak{D}}^{{\ell }_1}\left({}^C{\mathfrak{D}}^{{\ell }_2}\mu \right)\left(\mathfrak{t}\right)=K\left(\mathfrak{t},\nu \left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(\mathfrak{t}\right)\right)+\phi \left(\mathfrak{t}\right),(\mathfrak{t}\in \mathcal{O}),\hfill \\ {}^C{\mathfrak{D}}^{{\ell }_1^{*}}\left({}^C{\mathfrak{D}}^{{\ell }_2^{*}}\nu \right)\left(\mathfrak{t}\right)=M\left(\mathfrak{t},\mu \left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(\mathfrak{t}\right)\right)+\chi \left(\mathfrak{t}\right),(\mathfrak{t}\in \mathcal{O}),\hfill \\ \gamma \mu \left(0\right)=\delta \mu \left(1\right)=\alpha {}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(0\right)=\beta {}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(1\right)=0,\hfill \\ {\gamma }^{*}\nu \left(0\right)={\delta }^{*}\nu \left(1\right)={\alpha }^{*}{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\nu \left(0\right)={\beta }^{*}{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\nu \left(1\right)=0.\hfill \end{array}\right.$$

(47)

By Proposition (3), the solution of Equation (47) is given by the following.

$$\left\{\begin{array}{c}{\displaystyle \mu \left(\mathfrak{t}\right)={\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)\mathrm{d}r-\mathfrak{t}{\int }_0^1{\displaystyle \frac{{(1-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)\mathrm{d}r}\hfill \\ {\displaystyle +{\displaystyle \frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+2)}}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}K\left(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right)\right)\mathrm{d}r+{\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}\phi \left(r\right)\mathrm{d}r,}\hfill \\ {\displaystyle -\mathfrak{t}{\int }_0^1{\displaystyle \frac{{(1-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}\phi \left(r\right)\mathrm{d}r+{\displaystyle \frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+2)}}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}\phi \left(r\right)\mathrm{d}r,}\hfill \\ {\displaystyle \nu \left(\mathfrak{t}\right)={\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1^{*}+{\ell }_2^{*}-1}}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*})}}M\left(r,\mu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(r\right)\right)\mathrm{d}r-\mathfrak{t}{\int }_0^1{\displaystyle \frac{{(1-r)}^{{\ell }_1^{*}+{\ell }_2^{*}-1}}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*})}}M\left(r,\mu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(r\right)\right)\mathrm{d}r}\hfill \\ {\displaystyle +{\displaystyle \frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2^{*}+1}}{\mathsf{\Gamma }({\ell }_2^{*}+2)}}{\int }_0^1\frac{{(1-r)}^{{\ell }_1^{*}-1}}{\mathsf{\Gamma }\left({\ell }_1^{*}\right)}M\left(r,\mu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(r\right)\right)\mathrm{d}r+{\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1^{*}+{\ell }_2^{*}-1}}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*})}}\chi \left(r\right)\mathrm{d}r,}\hfill \\ {\displaystyle -\mathfrak{t}{\int }_0^1{\displaystyle \frac{{(1-r)}^{{\ell }_1^{*}+{\ell }_2^{*}-1}}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*})}}\chi \left(r\right)\mathrm{d}r+{\displaystyle \frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2^{*}+1}}{\mathsf{\Gamma }({\ell }_2^{*}+2)}}{\int }_0^1\frac{{(1-r)}^{{\ell }_1^{*}-1}}{\mathsf{\Gamma }\left({\ell }_1^{*}\right)}\chi \left(r\right)\mathrm{d}r.}\hfill \end{array}\right.$$

Since $\mathfrak{t}\in \mathcal{O},$ then from the last system, we have the following:

$$\begin{array}{cc}\hfill \left|\mu \right(\mathfrak{t})& {\displaystyle -[{\int }_0^{\mathfrak{t}}\frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))\mathrm{d}r-\mathfrak{t}{\int }_0^1\frac{{(1-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}K(r,\nu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))\mathrm{d}r}\hfill \\ \hfill & {\displaystyle +\frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}K(r,\mu \left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(r\right))\mathrm{d}r\left]\right|}\hfill \\ \hfill =& {\displaystyle |{\int }_0^{\mathfrak{t}}\frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}\phi \left(r\right)\mathrm{d}r-\mathfrak{t}{\int }_0^1\frac{{(1-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}\phi \left(r\right)\mathrm{d}r+\frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}\phi \left(r\right)\mathrm{d}r|}\hfill \\ \hfill \le & {\displaystyle {\int }_0^{\mathfrak{t}}\frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}|\phi \left(r\right)|\mathrm{d}r-\mathfrak{t}{\int }_0^1\frac{{(1-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}|\phi \left(r\right)|\mathrm{d}r+\frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+2)}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}|\phi \left(r\right)|\mathrm{d}r}\hfill \\ \hfill \le & \left(\frac{2}{\mathsf{\Gamma }({\ell }_1+{\ell }_2+1)}+\frac{1}{\mathsf{\Gamma }({\ell }_2+2)\mathsf{\Gamma }({\ell }_1+1)}\right){ϵ}_1={\eta }_1{ϵ}_1,\mathfrak{t}\in \mathcal{O},\hfill \end{array}$$

which provides the following.

$$|\mu \left(\mathfrak{t}\right)-H_1\nu \left(\mathfrak{t}\right)|\le {\eta }_1{ϵ}_1,\mathfrak{t}\in \mathcal{O}.$$

(48)

With the same computation techniques, we obtain the following.

$$|\nu \left(\mathfrak{t}\right)-H_2\mu \left(\mathfrak{t}\right)|\le {\eta }_2{ϵ}_2,\mathfrak{t}\in \mathcal{O}.$$

(49)

On the other hand, we have the following.

$$\begin{array}{c}\hfill |{}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(\mathfrak{t}\right)-{}^C{\mathfrak{D}}^{{\ell }_2}{H}_1\nu \left(\mathfrak{t}\right)|\le \frac{1}{\mathsf{\Gamma }(2-{\ell }_2)}{\int }_0^{\mathfrak{t}}{(\mathfrak{t}-r)}^{1-{\ell }_2}|{\mu }^{″}\left(r\right)-{\left(H_1\nu \right)}^{″}\left(r\right)|\mathrm{d}r\le \frac{{\gamma }_1}{\mathsf{\Gamma }(3-{\ell }_2)},\end{array}$$

(50)

With the same arguments, we obtain the following.

$$\begin{array}{cc}\hfill |{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\nu \left(\mathfrak{t}\right)-{}^C{\mathfrak{D}}^{{\ell }_2^{*}}{H}_2\mu \left(\mathfrak{t}\right)|\le & \frac{{\gamma }_2}{\mathsf{\Gamma }(3-{\ell }_2^{*})}.\hfill \end{array}$$

(51)

Therefore, from Equations (48)–(51), it follows that the following is the case.

$$\left\{\begin{array}{c}\parallel \mu -H_1\nu {\parallel }_{\mathfrak{B}}\le {\eta }_1{ϵ}_1+{\displaystyle \frac{{\gamma }_1}{\mathsf{\Gamma }(3-{\ell }_2)}},\hfill \\ \parallel \nu -H_2\mu {\parallel }_{\mathfrak{B}}\le {\eta }_2{ϵ}_2+{\displaystyle \frac{{\gamma }_2}{\mathsf{\Gamma }(3-{\ell }_2^{*})}},\hfill \end{array}\right.$$

The proof is completed. □

Theorem 5.

Consider the assumptions (H2) and (H3). If the following is the case:

$$\begin{array}{ccc}& & max\{\left[{\eta }_1{ϵ}_1+{\mathbf{C}}_K{\mathsf{\Delta }}_1{\eta }_2{ϵ}_2+\frac{{\gamma }_1}{\mathsf{\Gamma }(3-{\ell }_2)}+\frac{{\gamma }_2{\mathbf{C}}_K{\mathsf{\Delta }}_1}{\mathsf{\Gamma }(3-{\ell }_2^{*})}\right]{\left(1-{\mathbf{C}}_K{\mathbf{C}}_M{\mathsf{\Delta }}_1{\mathsf{\Delta }}_2\right)}^{-1},\hfill \\ & & \left[{\eta }_2{ϵ}_2+{\mathbf{C}}_M{\mathsf{\Delta }}_2{\eta }_1{ϵ}_1+\frac{{\gamma }_2}{\mathsf{\Gamma }(3-{\ell }_2^{*})}+\frac{{\gamma }_1{\mathbf{C}}_M{\mathsf{\Delta }}_2}{\mathsf{\Gamma }(3-{\ell }_2)}\right]{\left(1-{\mathbf{C}}_K{\mathbf{C}}_M{\mathsf{\Delta }}_1{\mathsf{\Delta }}_2\right)}^{-1}\}<1,\hfill \end{array}$$

where ${\mathbf{C}}_K{\mathbf{C}}_M{\mathsf{\Delta }}_1{\mathsf{\Delta }}_2<1,$ then the solution of the coupled BVPs of fractional Navier DEs (9) is H-U-stable.

Proof. 

Let $(\mu ,\nu )\in {\left(C(\mathcal{O},\mathbb{R})\right)}^2$ be a solution of the following system of inequalities:

$$\left\{\begin{array}{c}|{}^C{\mathfrak{D}}^{{\ell }_1}\left({}^C{\mathfrak{D}}^{{\ell }_2}\mu \right)\left(\mathfrak{t}\right)-K(\mathfrak{t},\nu \left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2}\nu \left(\mathfrak{t}\right))|\le {ϵ}_1,\mathfrak{t}\in \mathcal{O},\hfill \\ |{}^C{\mathfrak{D}}^{{\ell }_1^{*}}\left({}^C{\mathfrak{D}}^{{\ell }_2^{*}}\nu \right)\left(\mathfrak{t}\right)-M(\mathfrak{t},\mu \left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\mu \left(\mathfrak{t}\right))|\le {ϵ}_2,\mathfrak{t}\in \mathcal{O},\hfill \end{array}\right.$$

and $(\widehat{h},\tilde{h})\in {\left(C(\mathcal{O},\mathbb{R})\right)}^2$ be a unique solution of the following.

$$\left\{\begin{array}{c}{}^C{\mathfrak{D}}^{{\ell }_1}\left({}^C{\mathfrak{D}}^{{\ell }_2}\widehat{h}\right)\left(\mathfrak{t}\right)=K\left(\mathfrak{t},\tilde{h}\left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2}\tilde{h}\left(\mathfrak{t}\right)\right),\mathfrak{t}\in \mathcal{O},\hfill \\ {}^C{\mathfrak{D}}^{{\ell }_1^{*}}\left({}^C{\mathfrak{D}}^{{\ell }_2^{*}}\tilde{h}\right)\left(\mathfrak{t}\right)=M\left(\mathfrak{t},\widehat{h}\left(\mathfrak{t}\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\widehat{h}\left(\mathfrak{t}\right)\right),\mathfrak{t}\in \mathcal{O},\hfill \\ \gamma \widehat{h}\left(0\right)=\delta \widehat{h}\left(1\right)=\alpha {}^C{\mathfrak{D}}^{{\ell }_2}\widehat{h}\left(0\right)=\beta {}^C{\mathfrak{D}}^{{\ell }_2}\widehat{h}\left(1\right)=0,\hfill \\ {\gamma }^{*}\tilde{h}\left(0\right)={\delta }^{*}\tilde{h}\left(1\right)={\alpha }^{*}{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\tilde{h}\left(0\right)={\beta }^{*}{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\tilde{h}\left(1\right)=0.\hfill \end{array}\right.$$

(52)

Then, by Proposition (3) together with Equation (9), the solution of Equation (52) is given by

$$\begin{array}{cc}\hfill \widehat{h}\left(\mathfrak{t}\right)& {\displaystyle ={\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}K\left(r,\tilde{h}\left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\tilde{h}\left(r\right)\right)\mathrm{d}r-\mathfrak{t}{\int }_0^1{\displaystyle \frac{{(1-r)}^{{\ell }_1+{\ell }_2-1}}{\mathsf{\Gamma }({\ell }_1+{\ell }_2)}}K\left(r,\tilde{h}\left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\tilde{h}\left(r\right)\right)\mathrm{d}r}\hfill \\ \hfill & {\displaystyle +{\displaystyle \frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2+1}}{\mathsf{\Gamma }({\ell }_2+2)}}{\int }_0^1\frac{{(1-r)}^{{\ell }_1-1}}{\mathsf{\Gamma }\left({\ell }_1\right)}K\left(r,\tilde{h}\left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2}\tilde{h}\left(r\right)\right)\mathrm{d}r=H_1\tilde{h}\left(\mathfrak{t}\right),}\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill \tilde{h}\left(\mathfrak{t}\right)& {\displaystyle ={\int }_0^{\mathfrak{t}}{\displaystyle \frac{{(\mathfrak{t}-r)}^{{\ell }_1^{*}+{\ell }_2^{*}-1}}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*})}}M\left(r,\widehat{h}\left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\widehat{h}\left(r\right)\right)\mathrm{d}r-\mathfrak{t}{\int }_0^1{\displaystyle \frac{{(1-r)}^{{\ell }_1^{*}+{\ell }_2^{*}-1}}{\mathsf{\Gamma }({\ell }_1^{*}+{\ell }_2^{*})}}M\left(r,\widehat{h}\left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\widehat{h}\left(r\right)\right)\mathrm{d}r}\hfill \\ \hfill & {\displaystyle +{\displaystyle \frac{\mathfrak{t}-{\mathfrak{t}}^{{\ell }_2^{*}+1}}{\mathsf{\Gamma }({\ell }_2^{*}+2)}}{\int }_0^1\frac{{(1-r)}^{{\ell }_1^{*}-1}}{\mathsf{\Gamma }\left({\ell }_1^{*}\right)}M\left(r,\widehat{h}\left(r\right),{}^C{\mathfrak{D}}^{{\ell }_2^{*}}\widehat{h}\left(r\right)\right)\mathrm{d}r=H_2\widehat{h}\left(\mathfrak{t}\right).}\hfill \end{array}$$

Then, from the last system, we can write the following.

$$\begin{array}{ccc}\hfill |\mu \left(\mathfrak{t}\right)-\widehat{h}\left(\mathfrak{t}\right)|+|{}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(\mathfrak{t}\right)-{}^C{\mathfrak{D}}^{{\ell }_2}\widehat{h}\left(\mathfrak{t}\right)|& =& |\mu \left(\mathfrak{t}\right)-H_1\tilde{h}\left(\mathfrak{t}\right)|+|{}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(\mathfrak{t}\right)-{}^C{\mathfrak{D}}^{{\ell }_2}{H}_1\tilde{h}\left(\mathfrak{t}\right)|\hfill \\ & =& |\mu \left(\mathfrak{t}\right)-H_1\nu \left(\mathfrak{t}\right)+H_1\nu \left(\mathfrak{t}\right)-H_1\tilde{h}\left(\mathfrak{t}\right)|\hfill \\ \hfill \\ & +& |{}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(\mathfrak{t}\right)-{}^C{\mathfrak{D}}^{{\ell }_2}{H}_1\nu \left(\mathfrak{t}\right)+{}^C{\mathfrak{D}}^{{\ell }_2}{H}_1\nu \left(\mathfrak{t}\right)-{}^C{\mathfrak{D}}^{{\ell }_2}{H}_1\tilde{h}\left(\mathfrak{t}\right)|\hfill \\ \hfill \\ & \le & |\mu \left(\mathfrak{t}\right)-H_1\nu \left(\mathfrak{t}\right)|+|H_1\nu \left(\mathfrak{t}\right)-H_1\tilde{h}\left(\mathfrak{t}\right)|\hfill \\ \hfill \\ & +& |{}^C{\mathfrak{D}}^{{\ell }_2}\mu \left(\mathfrak{t}\right)-{}^C{\mathfrak{D}}^{{\ell }_2}{H}_1\nu \left(\mathfrak{t}\right)|+|{}^C{\mathfrak{D}}^{{\ell }_2}{H}_1\nu \left(\mathfrak{t}\right)-{}^C{\mathfrak{D}}^{{\ell }_2}{H}_1\tilde{h}\left(\mathfrak{t}\right)|.\hfill \end{array}$$

From the last inequalities and Lemma (1) with some modifications, we obtain the following.

$$\begin{array}{ccc}\hfill \parallel \mu -\widehat{h}{\parallel }_{\mathfrak{B}}& \le & \parallel \mu -H_1\nu {\parallel }_{\mathfrak{B}}+\parallel H_1\nu -H_1\tilde{h}{\parallel }_{\mathfrak{B}}\hfill \\ & \le & {\eta }_1{\epsilon }_1+\frac{{\gamma }_1}{\mathsf{\Gamma }(3-{\ell }_2)}+{\mathbf{C}}_K{\mathsf{\Delta }}_1\parallel \nu -\tilde{h}{\parallel }_{\mathfrak{B}}.\hfill \end{array}$$

(53)

In an analogous manner, we find the following.

$$\begin{array}{c}\hfill \parallel \nu -\tilde{h}{\parallel }_{\mathfrak{B}}\le {\eta }_2{\epsilon }_2+\frac{{\gamma }_2}{\mathsf{\Gamma }(3-{\ell }_2^{*})}+{\mathbf{C}}_M{\mathsf{\Delta }}_2\parallel \mu -\widehat{h}{\parallel }_{\mathfrak{B}}.\end{array}$$

(54)

Since we have assumed that ${\mathbf{C}}_K{\mathbf{C}}_M{\mathsf{\Delta }}_1{\mathsf{\Delta }}_2<1,$ then from Equations (53) and (54), we obtain

$$\parallel \mu -\widehat{h}{\parallel }_{\mathfrak{B}} \le \left[{\eta }_1{ϵ}_1+{\mathbf{C}}_K{\mathsf{\Delta }}_1{\eta }_2{ϵ}_2+\frac{{\gamma }_1}{\mathsf{\Gamma }(3-{\ell }_2)}+\frac{{\gamma }_2{\mathbf{C}}_K{\mathsf{\Delta }}_1}{\mathsf{\Gamma }(3-{\ell }_2^{*})}\right]{\left(1-{\mathbf{C}}_K{\mathbf{C}}_M{\mathsf{\Delta }}_1{\mathsf{\Delta }}_2\right)}^{-1},$$

(55)

and

$$\parallel \nu -\tilde{h}{\parallel }_{\mathfrak{B}} \le \left[{\eta }_2{ϵ}_2+{\mathbf{C}}_M{\mathsf{\Delta }}_2{\eta }_1{ϵ}_1+\frac{{\gamma }_2}{\mathsf{\Gamma }(3-{\ell }_2^{*})}+\frac{{\gamma }_1{\mathbf{C}}_M{\mathsf{\Delta }}_2}{\mathsf{\Gamma }(3-{\ell }_2)}\right]{\left(1-{\mathbf{C}}_K{\mathbf{C}}_M{\mathsf{\Delta }}_1{\mathsf{\Delta }}_2\right)}^{-1}.$$

(56)

Consequently, Equations (55) and (56) result in the following.

$$\begin{array}{ccc}\hfill \parallel (\mu ,\nu )-\left(\widehat{h},\tilde{h}\right){\parallel }_{\mathfrak{B}\times \mathfrak{B}}& \le & max\{\left[{\eta }_1{ϵ}_1+{\mathbf{C}}_K{\mathsf{\Delta }}_1{\eta }_2{ϵ}_2+\frac{{\gamma }_1}{\mathsf{\Gamma }(3-{\ell }_2)}+\frac{{\gamma }_2{\mathbf{C}}_K{\mathsf{\Delta }}_1}{\mathsf{\Gamma }(3-{\ell }_2^{*})}\right]{\left(1-{\mathbf{C}}_K{\mathbf{C}}_M{\mathsf{\Delta }}_1{\mathsf{\Delta }}_2\right)}^{-1},\hfill \\ & & \left[{\eta }_2{ϵ}_2+{\mathbf{C}}_M{\mathsf{\Delta }}_2{\eta }_1{ϵ}_1+\frac{{\gamma }_2}{\mathsf{\Gamma }(3-{\ell }_2^{*})}+\frac{{\gamma }_1{\mathbf{C}}_M{\mathsf{\Delta }}_2}{\mathsf{\Gamma }(3-{\ell }_2)}\right]{\left(1-{\mathbf{C}}_K{\mathbf{C}}_M{\mathsf{\Delta }}_1{\mathsf{\Delta }}_2\right)}^{-1}\}.\hfill \end{array}$$

Consequently, the solution of the coupled BVPs of fractional Navier differential Equations (9) is H-U-stable. □

5. Numerical Solutions via GDT-Method

Due to the complexity of the nonlinear FDEs, it would be time consuming to obtain the exact solutions of given boundary problems with complicated nonlinear boundary conditions. Hence, we need to try deriving new approximate techniques and methods to solve these nonlinear problems. In the existing techniques, there exist some numerical methods that are applicable in finding the approximate solutions for nonlinear BVPs of non-integer orders. In the present study, we aim to apply one of these methods called the differential transformation method, which Zhou has presented in [38] and Odibat et al. had later extended it to the generalized version called GDT-Method in [39]. It is an iterative technique that yields analytical solutions of the Taylor series to FDEs with arbitrary boundary or initial conditions. In fact, it is with respect to the generalized Taylor’s formulas. For this reason, we shall use the GDT-Method to obtain approximate solutions of our coupled Navier FBVPs (1).

The generalized differential transform for the ${\mathfrak{s}}^{th}$-derivative of $\mu \left(\mathfrak{t}\right)$ and $\nu \left(\mathfrak{t}\right)$ in terms of one variable is introduced by

$${\mu }^{🟉}\left(\mathfrak{s}\right)={\displaystyle \frac{1}{\mathsf{\Gamma }(1+\mathfrak{s}\tau )}}\left[{\left({}^c{\mathfrak{D}}^{\tau }\right)}^{\mathfrak{s}}\mu \left(\mathfrak{t}\right)\right]{|}_{\mathfrak{t}=0},{\nu }^{🟉}\left(\mathfrak{s}\right)={\displaystyle \frac{1}{\mathsf{\Gamma }(1+\mathfrak{s}\tau )}}\left[{\left({}^c{\mathfrak{D}}^{\tau }\right)}^{\mathfrak{s}}\nu \left(\mathfrak{t}\right)\right]{|}_{\mathfrak{t}=0},$$

(57)

in which ${\left({}^c{\mathfrak{D}}^{\tau }\right)}^{\mathfrak{s}}=\stackrel{\mathfrak{s}-times}{\overbrace{{}^c{\mathfrak{D}}^{\tau }.{}^c{\mathfrak{D}}^{\tau }\cdots {}^c{\mathfrak{D}}^{\tau }}}$ [40]. The inverses are given as follows.

$${\displaystyle \mu \left(\mathfrak{t}\right)=\sum _{\mathfrak{s}=0}^{\infty }{\mu }^{🟉}\left(\mathfrak{s}\right){\mathfrak{t}}^{\tau \mathfrak{s}},\nu \left(\mathfrak{t}\right)=\sum _{\mathfrak{s}=0}^{\infty }{\nu }^{🟉}\left(\mathfrak{s}\right){\mathfrak{t}}^{\tau \mathfrak{s}}.}$$

Therefore, the approximate solutions of the Navier model (1) are found in the form of a finite series of the analytical polynomials:

$${\displaystyle \mu \left(\mathfrak{t}\right)=\sum _{\mathfrak{s}=0}^N{\mu }^{🟉}\left(\mathfrak{s}\right){\mathfrak{t}}^{\tau \mathfrak{s}},\nu \left(\mathfrak{t}\right)=\sum _{\mathfrak{s}=0}^N{\nu }^{🟉}\left(\mathfrak{s}\right){\mathfrak{t}}^{\tau \mathfrak{s}},}$$

(58)

where $\tau$ stands for the order of the mentioned transformation and must be taken so that $\mathfrak{l}\tau ={\ell }_2,$ $\mathfrak{p}\tau ={\ell }_2^{*},$ $\mathfrak{m}\tau ={\ell }_1+{\ell }_2$, and $\mathfrak{n}\tau ={\ell }_1^{*}+{\ell }_2^{*}$ with $\mathfrak{l},\mathfrak{p},\mathfrak{m},\mathfrak{n}\in \mathbb{N}$ and ${\mu }^{🟉}\left(\mathfrak{s}\right),$${\nu }^{🟉}\left(\mathfrak{s}\right)$ are the generalized differential transforms of $\mu \left(\mathfrak{t}\right)$ and $\nu \left(\mathfrak{t}\right)$, respectively, and they are given as

$$\left\{\begin{array}{c}{\displaystyle {\mu }^{🟉}(\mathfrak{m}+\mathfrak{s})=\frac{\mathsf{\Gamma }((\mathfrak{m}+\mathfrak{s})\tau -{\ell }_1-{\ell }_2+1)}{\mathsf{\Gamma }\left(\right(\mathfrak{m}+\mathfrak{s})\tau +1)}{K}^{🟉}(\mathfrak{s},{\nu }^{🟉}\left(\mathfrak{s}\right)),}\hfill \\ {\displaystyle {\nu }^{🟉}(\mathfrak{n}+\mathfrak{s})=\frac{\mathsf{\Gamma }((\mathfrak{n}+\mathfrak{s})\tau -{\ell }_1^{*}-{\ell }_2^{*}+1)}{\mathsf{\Gamma }\left(\right(\mathfrak{n}+\mathfrak{s})\tau +1)}{M}^{🟉}(\mathfrak{s},{\mu }^{🟉}\left(\mathfrak{s}\right)),}\hfill \end{array}\right.$$

(59)

in which $K^{🟉}(\mathfrak{s},{\nu }^{🟉}\left(\mathfrak{s}\right))$ and $M^{🟉}(\mathfrak{s},{\mu }^{🟉}\left(\mathfrak{s}\right))$ are the generalized ${\tau }^{th}$-differential transforms of $K(\mathfrak{s},\nu (\mathfrak{s}\left)\right)$ and $M(s,\mu (\mathfrak{s}\left)\right)$, respectively.

Since $\mu \left(0\right)=0,$ $\nu \left(0\right)=0,$ ${}^c{\mathfrak{D}}^{{\ell }_2}\mu \left(0\right)=0$, and ${}^c{\mathfrak{D}}^{{\ell }_2^{*}}\nu \left(0\right)=0,$ so their GDTs provide ${\mu }^{🟉}\left(\mathfrak{l}\right)=\widehat{c},$ ${\nu }^{🟉}\left(\mathfrak{p}\right)=\tilde{c}$, ${\mu }^{🟉}\left(0\right)=0$, and ${\nu }^{🟉}\left(0\right)=0$, ${\mu }^{🟉}\left(\mathfrak{s}\right)=0,$ for all $\mathfrak{s}$ satisfying $0<\mathfrak{s}\tau <{\ell }_2$ or $\mathfrak{l}<\mathfrak{s}<\mathfrak{m}$, ${\nu }^{🟉}\left(\mathfrak{s}\right)=0,$ for all $\mathfrak{s}$ satisfying $0<\mathfrak{s}\tau <{\ell }_2^{*}$ or $\mathfrak{p}<\mathfrak{s}<\mathfrak{n}$, where $\widehat{c}$ and $\tilde{c}\in \mathbb{R}$ are unknowns that can be computed by applying other initial conditions of the Navier model (1).

By recursive Equation (59), the solution $\left(\mu \right(\mathfrak{t}),\nu (\mathfrak{t}\left)\right)$ of Navier problem (1) can be represented via the following series:

$${\displaystyle \mu \left(\mathfrak{t}\right)=\sum _{\mathfrak{s}=0}^N{\mu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right){\mathfrak{t}}^{\tau \mathfrak{s}},\nu \left(\mathfrak{t}\right)=\sum _{\mathfrak{s}=0}^N{\nu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right){\mathfrak{t}}^{\tau \mathfrak{s}},}$$

(60)

where ${\mu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right)$, and ${\nu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right)$ are coefficients in terms of $\widehat{c}$ and $\tilde{c}$, which can be obtained by invoking other initial conditions of given Navier problem (1).

From (60), we have the following.

$$\left\{\begin{array}{c}{\displaystyle \mu \left(1\right)=\sum _{\mathfrak{s}=0}^N{\mu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right),}\hfill \\ {\displaystyle \nu \left(1\right)=\sum _{\mathfrak{s}=0}^N{\nu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right),}\hfill \\ {\displaystyle \delta {}^c{\mathfrak{D}}^{{\ell }_2}\mu \left(1\right)=\frac{\beta }{\mathsf{\Gamma }(2-{\ell }_2)}{\int }_0^1{(1-\lambda )}^{1-{\ell }_2}\sum _{\mathfrak{s}=0}^N\mathfrak{s}{\mu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right)\tau s\left(\tau s-1\right){\lambda }^{\tau \mathfrak{s}-2}d\lambda ,}\hfill \\ {\displaystyle {\delta }^{*}{}^c{\mathfrak{D}}^{{\ell }_2}\nu \left(1\right)=\frac{{\beta }^{*}}{\mathsf{\Gamma }(2-{\ell }_2^{*})}{\int }_0^1{(1-\lambda )}^{1-{\ell }_2^{*}}\sum _{\mathfrak{s}=0}^N\mathfrak{s}{\nu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right)\tau s\left(\tau s-1\right){\lambda }^{\tau \mathfrak{s}-2}d\lambda .}\hfill \end{array}\right.$$

(61)

Therefore, Equation (61) yields the following system.

$$\left\{\begin{array}{c}{\displaystyle \delta \sum _{\mathfrak{s}=0}^N{\mu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right)-\frac{\beta {\tau }^2}{\mathsf{\Gamma }(2-{\ell }_2)}\sum _{\mathfrak{s}=0}^N{\mathfrak{s}}^2{\mu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right){\int }_0^1{(1-\lambda )}^{1-{\ell }_2}{\lambda }^{\tau \mathfrak{s}-2}d\lambda }\hfill \\ {\displaystyle +\frac{\beta \tau }{\mathsf{\Gamma }(2-{\ell }_2)}\sum _{\mathfrak{s}=0}^N\mathfrak{s}{\mu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right){\int }_0^1{(1-\lambda )}^{1-{\ell }_2}{\lambda }^{\tau \mathfrak{s}-2}d\lambda =0}\hfill \\ {\displaystyle {\delta }^{*}\sum _{\mathfrak{s}=0}^N{\nu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right)-\frac{{\beta }^{*}{\tau }^2}{\mathsf{\Gamma }(2-{\ell }_2^{*})}\sum _{\mathfrak{s}=0}^N{\mathfrak{s}}^2{\mu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right){\int }_0^1{(1-\lambda )}^{1-{\ell }_2^{*}}{\lambda }^{\tau \mathfrak{s}-2}d\lambda }\hfill \\ {\displaystyle +\frac{{\beta }^{*}\tau }{\mathsf{\Gamma }(2-{\ell }_2^{*})}\sum _{\mathfrak{s}=0}^N\mathfrak{s}{\nu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right){\int }_0^1{(1-\lambda )}^{1-{\ell }_2^{*}}{\lambda }^{\tau \mathfrak{s}-2}d\lambda =0}\hfill \end{array}\right.$$

(62)

Finally, by solving Equation (62) in terms of $\widehat{c}$ and $\tilde{c}$ and replacing them in Equation (60), we find the approximate solutions of Navier model (1).

Estimation of Error

The base of the GDT-method consists in finding a development in fractional series for solutions of nonlinear models containing fractional derivatives around the initial value ${\mathfrak{t}}_0$.

$$\left\{\begin{array}{c}{\displaystyle \mu \left(\mathfrak{t}\right)=\sum _{k=0}^{+\infty }{a}_k{(\mathfrak{t}-{\mathfrak{t}}_0)}^{\mathfrak{s}k},k\in I=({\mathfrak{t}}_0,{\mathfrak{t}}_0+r),r>0,}\hfill \\ {\displaystyle \nu \left(\mathfrak{t}\right)=\sum _{k=0}^{+\infty }{b}_k{(\mathfrak{t}-{\mathfrak{t}}_0)}^{\mathfrak{s}k},k\in I=({\mathfrak{t}}_0,{\mathfrak{t}}_0+r),r>0.}\hfill \end{array}\right.$$

(63)

For estimation of error, the following theorems are useful.

Theorem 6

([41]). Let ${\mathsf{\Psi }}_k\left(\mathfrak{t}\right)=a_k{(\mathfrak{t}-{\mathfrak{t}}_0)}^{\mathfrak{s}k}$. Then, the series solution ${\displaystyle \sum _{k=0}^{+\infty }{\mathsf{\Psi }}_k\left(\mathfrak{t}\right)}$ converges if there exists $0<\gamma <1$ such that the following is the case.

$$\parallel {\mathsf{\Psi }}_{k+1}\left(\mathfrak{t}\right)\parallel \le \gamma \parallel {\mathsf{\Psi }}_k\left(\mathfrak{t}\right)\parallel ,\forall k\ge k_0,for\hspace{1em}some k_0\in \mathbb{N}.$$

Theorem 7

([41]). Suppose that the series solution ${\displaystyle \sum _{k=0}^{+\infty }{\mathsf{\Psi }}_k\left(\mathfrak{t}\right)}$ with ${\mathsf{\Psi }}_k\left(\mathfrak{t}\right)=a_k{(\mathfrak{t}-{\mathfrak{t}}_0)}^{\mathfrak{s}k}$ converges to the solution $\mu \left(\mathfrak{t}\right)$. If the truncated series ${\displaystyle \sum _{k=0}^m{\mathsf{\Psi }}_k\left(\mathfrak{t}\right)}$ is an approximation to the solution $\mu \left(\mathfrak{t}\right)$, then the maximum absolute truncated error is estimated as follows.

$${\displaystyle \parallel \mu \left(\mathfrak{t}\right)-\sum _{k=0}^m{\mathsf{\Psi }}_k\left(\mathfrak{t}\right)\parallel \le \frac{1}{1-\gamma }{\gamma }^{m-m_0+1}\underset{\mathfrak{t}\in I}{max}|a_{m_0}{(\mathfrak{t}-{\mathfrak{t}}_0)}^{\mathfrak{s}{m}_0}|,for any{m}_0\ge 0and{a}_{m_0}\ne 0.}$$

In conclusion, in view of Theorem 6, we find that the fractional power series solution ${\displaystyle \mu \left(\mathfrak{t}\right)=\sum _{k=0}^{+\infty }{a}_k{(\mathfrak{t}-{\mathfrak{t}}_0)}^{\mathfrak{s}k}}$ converges to the exact solution if there exists $0<\gamma <1$ such that

$\parallel {\mathsf{\Psi }}_{k+1}\left(\mathfrak{t}\right)\parallel \le \gamma \parallel {\mathsf{\Psi }}_k\left(\mathfrak{t}\right)\parallel ,\forall k\ge k_0,\mathrm{for} \mathrm{some}{k}_0\in \mathbb{N}.$ In other words, if for each $j\ge k_0$, we have the folllowing:

$${\gamma }_{j+1}=\left\{\begin{array}{c}\frac{\parallel {\mathsf{\Psi }}_{j+1}\left(\mathfrak{t}\right)\parallel }{\parallel {\mathsf{\Psi }}_j\left(\mathfrak{t}\right)\parallel },\parallel {\mathsf{\Psi }}_j\left(\mathfrak{t}\right)\parallel \ne 0,\hfill \\ 0,\parallel {\mathsf{\Psi }}_j\left(\mathfrak{t}\right)\parallel =0,\hfill \end{array}\right.$$

$j\in \mathbb{N}\cup \left\{0\right\},$ where ${\displaystyle \parallel {\mathsf{\Psi }}_j\left(\mathfrak{t}\right)\parallel =\underset{\mathfrak{t}\in I}{max}|a_j{(\mathfrak{t}-{\mathfrak{t}}_0)}^{\mathfrak{s}j}|}$ then, the series solution ${\displaystyle \sum _{k=0}^{+\infty }{\mathsf{\Psi }}_k\left(\mathfrak{t}\right)}$ converges to the exact solution $\mu \left(\mathfrak{t}\right)$, when $0\le {\gamma }_j\le 1$, for all $j\ge k_0$. In addition, the conclusion of Theorem 7 confirms that the maximum absolute truncation error is estimated to be the following:

$${\displaystyle \parallel \mu \left(\mathfrak{t}\right)-\sum _{k=0}^m{\mathsf{\Psi }}_k\left(\mathfrak{t}\right)\parallel \le \frac{1}{1-\beta }{\gamma }^{m-m_0+1}\underset{\mathfrak{t}\in I}{max}|a_{m_0}{(\mathfrak{t}-{\mathfrak{t}}_0)}^{\mathfrak{s}{m}_0}|,}$$

where $\beta =\{{\gamma }_j,m_0+1,m_0+2,\dots ,m+1\}.$ Consequently, if the trucated series

$${\displaystyle \left(\sum _{k=0}^m{\mathsf{\Psi }}_k\left(\mathfrak{t}\right),\sum _{k=0}^n{\mathsf{\Phi }}_k\left(\mathfrak{t}\right)\right)}$$

is an approximation to the solution $\left(\mu \left(\mathfrak{t}\right),\nu \left(\mathfrak{t}\right)\right)$ defined by Equation (63), then the maximum absolute truncated error is estimated as follows:

$$\left\{\begin{array}{c}{\displaystyle \parallel \mu \left(\mathfrak{t}\right)-\sum _{k=0}^m{\mathsf{\Psi }}_k\left(\mathfrak{t}\right)\parallel \le \frac{1}{1-\beta }{\gamma }^{m-m_0+1}\underset{\mathfrak{t}\in I}{max}|a_{m_0}{(\mathfrak{t}-{\mathfrak{t}}_0)}^{\mathfrak{s}{m}_0}|,\mathrm{for} \mathrm{any}{m}_0\ge 0\mathrm{and}{a}_{m_0}\ne 0.}\hfill \\ {\displaystyle \parallel \nu \left(\mathfrak{t}\right)-\sum _{k=0}^n{\mathsf{\Phi }}_k\left(\mathfrak{t}\right)\parallel \le \frac{1}{1-\delta }{\delta }^{n-n_0+1}\underset{\mathfrak{t}\in I}{max}|a_{n_0}{(\mathfrak{t}-{\mathfrak{t}}_0)}^{\mathfrak{s}{n}_0}|,\mathrm{for} \mathrm{any}{n}_0\ge 0\mathrm{and}{b}_{n_0}\ne 0,}\hfill \end{array}\right.$$

where $\beta =\{{\gamma }_j,m_0+1,m_0+2,\dots ,m+1\}$ and $\delta =\{{\delta }_j,n_0+1,n_0+2,\dots ,n+1\}$.

6. Examples

Example 1.

We design the following system of coupled sequential Navier FBVPs:

$$\left\{\begin{array}{c}{\displaystyle {}^C{\mathfrak{D}}^{\frac{5}{4}}\left({}^C{\mathfrak{D}}^{\frac{7}{4}}\mu \right)\left(\mathfrak{t}\right)=\frac{e^{-\mathfrak{t}}\left(|\nu \left(\mathfrak{t}\right)|+|{}^C{\mathfrak{D}}^{\frac{7}{4}}\nu \left(\mathfrak{t}\right)|\right)}{\left(9+e^{-\mathfrak{t}}\right)\left(1+|\nu \left(\mathfrak{t}\right)|+|{}^C{\mathfrak{D}}^{\frac{7}{4}}\nu \left(\mathfrak{t}\right)|\right)},(\mathfrak{t}\in \mathcal{O}:=[0,1]),}\hfill \\ {\displaystyle {}^C{\mathfrak{D}}^{\frac{3}{2}}\left({}^C{\mathfrak{D}}^{\frac{5}{4}}\nu \right)\left(\mathfrak{t}\right)=\frac{|\mu \left(\mathfrak{t}\right)|+|{}^C{\mathfrak{D}}^{\frac{8}{5}}\mu \left(\mathfrak{t}\right)|}{\left(\pi +e^{-2\mathfrak{t}}\right)\left(1+|\mu \left(\mathfrak{t}\right)|+|{}^C{\mathfrak{D}}^{\frac{8}{5}}\mu \left(\mathfrak{t}\right)|\right)},(\mathfrak{t}\in \mathcal{O}:=[0,1]),}\hfill \\ \mu \left(0\right)=\mu \left(1\right)={}^C{\mathfrak{D}}^{\frac{7}{4}}\mu \left(0\right)={}^C{\mathfrak{D}}^{\frac{7}{4}}\mu \left(1\right)=0,\hfill \\ \nu \left(0\right)=\nu \left(1\right)={}^C{\mathfrak{D}}^{\frac{5}{4}}\nu \left(0\right)={}^C{\mathfrak{D}}^{\frac{5}{4}}\nu \left(1\right)=0,\hfill \end{array}\right.$$

(64)

where ${\ell }_1={\displaystyle \frac{5}{4}},$ ${\ell }_2={\displaystyle \frac{7}{4}},$ ${\ell }_1^{*}={\displaystyle \frac{3}{2}},$ ${\ell }_2^{*}={\displaystyle \frac{5}{4}}$, and the following is the case.

$${\displaystyle K(\mathfrak{t},{\nu }_1,{\nu }_2)=\frac{e^{-\mathfrak{t}}\left(|{\nu }_1\left(\mathfrak{t}\right)|+|{\nu }_2\left(\mathfrak{t}\right)|\right)}{\left(9+e^{-\mathfrak{t}}\right)\left(1+|{\nu }_1\left(\mathfrak{t}\right)|+|{\nu }_2\left(\mathfrak{t}\right)|\right)},M(\mathfrak{t},{\mu }_1,{\mu }_2)=\frac{|{\mu }_1\left(\mathfrak{t}\right)|+|{\mu }_2\left(\mathfrak{t}\right)|}{\left({\pi }^2+e^{-2\mathfrak{t}}\right)\left(1+|{\mu }_1\left(\mathfrak{t}\right)|+|{\mu }_2\left(\mathfrak{t}\right)|\right)}.}$$

Then, for each $\mathfrak{t}\in [0,1]$ and ${\mu }_1,{\mu }_2,{\widehat{\mu }}_1,{\widehat{\mu }}_2,{\nu }_1,{\nu }_2,{\widehat{\nu }}_1,{\widehat{\nu }}_2\in [0,+\infty ),$ we have

$$\begin{array}{ccc}\hfill |K(\mathfrak{t},{\nu }_1,{\nu }_2)-K(\mathfrak{t},{\widehat{\nu }}_1,{\widehat{\nu }}_2)|& =& \frac{e^{-\mathfrak{t}}}{9+e^{-\mathfrak{t}}}|\frac{{\nu }_1\left(\mathfrak{t}\right)+{\nu }_2\left(\mathfrak{t}\right)}{1+{\nu }_1\left(\mathfrak{t}\right)+{\nu }_2\left(\mathfrak{t}\right)}-\frac{{\widehat{\nu }}_1\left(\mathfrak{t}\right)+{\widehat{\nu }}_2\left(\mathfrak{t}\right)}{1+{\widehat{\nu }}_1\left(\mathfrak{t}\right)+{\widehat{\nu }}_2\left(\mathfrak{t}\right)}|\hfill \\ & =& \frac{e^{-\mathfrak{t}}}{9+e^{-\mathfrak{t}}}\left|\frac{{\nu }_1\left(\mathfrak{t}\right)-{\widehat{\nu }}_1\left(\mathfrak{t}\right)+{\nu }_2\left(\mathfrak{t}\right)-{\widehat{\nu }}_2\left(\mathfrak{t}\right)}{\left(1+{\nu }_1\left(\mathfrak{t}\right)+{\nu }_2\left(\mathfrak{t}\right)\right)\left(1+{\widehat{\nu }}_1\left(\mathfrak{t}\right)+{\widehat{\nu }}_2\left(\mathfrak{t}\right)\right)}\right|\hfill \\ & \le & \frac{1}{9}\left[|{\nu }_1\left(\mathfrak{t}\right)-{\widehat{\nu }}_1\left(\mathfrak{t}\right)|+|{\nu }_2\left(\mathfrak{t}\right)-{\widehat{\nu }}_2\left(\mathfrak{t}\right)|\right],\hfill \end{array}$$

and

$$\begin{array}{c}\hfill |M(\mathfrak{t},{\mu }_1,{\mu }_2)-M(\mathfrak{t},{\widehat{\mu }}_1,{\widehat{\mu }}_2)|\le \frac{1}{{\pi }^2}\left[|{\mu }_1\left(\mathfrak{t}\right)-{\widehat{\mu }}_1\left(\mathfrak{t}\right)|+|{\mu }_2\left(\mathfrak{t}\right)-{\widehat{\mu }}_2\left(\mathfrak{t}\right)|\right].\end{array}$$

Thus, ${\mathbf{C}}_K={\displaystyle \frac{1}{9}},$ ${\mathbf{C}}_M={\displaystyle \frac{1}{{\pi }^2}},$ ${\mathsf{\Delta }}_1\approx 2.6662$ and ${\mathsf{\Delta }}_2\approx 2.8342.$ Hence, ${\mathbf{C}}_K{\mathsf{\Delta }}_1<1,$ ${\mathbf{C}}_M{\mathsf{\Delta }}_2<1.$

We observe that Theorem (3) is valid. Consequently, the system of sequential Navier FBVPs (64) has one solution.

Here, we try to find an approximate solution of the given sequential Navier FBVP (64) via the GDT-Method. Choose $\tau ={\displaystyle \frac{1}{4}}$, which yields $\mathfrak{l}=7,$ $\mathfrak{p}=5,$ $\mathfrak{m}=12$ and $\mathfrak{n}=11.$ By virtue of Equation (59) for the sequential Navier FBVPs (64), we have the following.

$$\left\{\begin{array}{c}{\displaystyle {\mu }^{🟉}(\mathfrak{s}+12)=\frac{\mathsf{\Gamma }\left(\frac{\mathfrak{s}}{4}+1\right)}{\mathsf{\Gamma }\left(\frac{\mathfrak{s}}{4}+4\right)}{K}^{🟉}(\mathfrak{s},{\nu }^{🟉}\left(\mathfrak{s}\right)),}\hfill \\ {\displaystyle {\nu }^{🟉}(\mathfrak{s}+11)=\frac{\mathsf{\Gamma }\left(\frac{\mathfrak{s}}{4}+1\right)}{\mathsf{\Gamma }\left(\frac{\mathfrak{s}}{4}+\frac{15}{4}\right)}{M}^{🟉}(\mathfrak{s},{\mu }^{🟉}\left(\mathfrak{s}\right)).}\hfill \end{array}\right.$$

(65)

Then, we obtain the following.

$$\begin{array}{cccccc}{\mu }^{🟉}\left(0\right)=0,& {\mu }^{🟉}\left(1\right)=0,& {\mu }^{🟉}\left(2\right)=0,& {\mu }^{🟉}\left(3\right)=0,& {\mu }^{🟉}\left(4\right)=0,& {\mu }^{🟉}\left(5\right)=0,\\ {\mu }^{🟉}\left(6\right)=0,& {\mu }^{🟉}\left(7\right)=\widehat{c},& {\mu }^{🟉}\left(8\right)=0,& {\mu }^{🟉}\left(9\right)=0,& {\mu }^{🟉}\left(10\right)=0,& {\mu }^{🟉}\left(11\right)=0,\\ {\nu }^{🟉}\left(0\right)=0,& {\nu }^{🟉}\left(1\right)=0,& {\nu }^{🟉}\left(2\right)=0,& {\nu }^{🟉}\left(3\right)=0,& {\nu }^{🟉}\left(4\right)=0,& {\nu }^{🟉}\left(5\right)=\tilde{c},\\ {\nu }^{🟉}\left(6\right)=0,& {\nu }^{🟉}\left(7\right)=\widehat{c},& {\nu }^{🟉}\left(8\right)=0,& {\nu }^{🟉}\left(9\right)=0,& {\nu }^{🟉}\left(10\right)=0.\end{array}$$

By Equation (62), we obtain the following.

$$\left\{\begin{array}{cc}{\displaystyle \sum _{\mathfrak{s}=0}^N{\mu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right)}\hfill & {\displaystyle -\frac{1}{16\mathsf{\Gamma }\left(0.25\right)}\sum _{\mathfrak{s}=0}^N{\mathfrak{s}}^2{\mu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right){\int }_0^1{(1-\lambda )}^{-0.75}{\lambda }^{\frac{s}{4}-2}d\lambda }\hfill \\ \hfill & {\displaystyle +\frac{1}{4\mathsf{\Gamma }\left(0.25\right)}\sum _{\mathfrak{s}=0}^N\mathfrak{s}{\mu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right){\int }_0^1{(1-\lambda )}^{-0.75}{\lambda }^{\frac{s}{4}-2}d\lambda =0}\hfill \\ {\displaystyle \sum _{\mathfrak{s}=0}^N{\nu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right)}\hfill & {\displaystyle -\frac{1}{16\mathsf{\Gamma }\left(0.75\right)}\sum _{\mathfrak{s}=0}^N{\mathfrak{s}}^2{\mu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right){\int }_0^1{(1-\lambda )}^{-0.25}{\lambda }^{\frac{s}{4}-2}d\lambda }\hfill \\ \hfill & {\displaystyle +\frac{1}{4\mathsf{\Gamma }\left(0.75\right)}\sum _{\mathfrak{s}=0}^N\mathfrak{s}{\nu }_{\widehat{c},\tilde{c}}^{🟉}\left(\mathfrak{s}\right){\int }_0^1{(1-\lambda )}^{-0.25}{\lambda }^{\frac{s}{4}-2}d\lambda =0.}\hfill \end{array}\right.$$

(66)

By using the recurrence relationship (65) truncated, respectively, at $\mathfrak{s}=15$ and $\mathfrak{s}=20$ and computing the constants $\widehat{c}$ and $\tilde{c}$ from Equation (66), the approximate solutions $({\mu }_1\left(\mathfrak{t}\right),{\nu }_1\left(\mathfrak{t}\right)),$ $({\mu }_2\left(\mathfrak{t}\right)$, and ${\nu }_2\left(\mathfrak{t}\right))$ are obtained by

$$\left\{\begin{array}{c}{\mu }_1\left(\mathfrak{t}\right)=0.1103{\mathfrak{t}}^{1.75}-0.1247{\mathfrak{t}}^3+0.1054{\mathfrak{t}}^{3.25}-0.0063{\mathfrak{t}}^{3.5}+0.0054{\mathfrak{t}}^{3.75},\hfill \\ {\nu }_1\left(\mathfrak{t}\right)=0.7307{\mathfrak{t}}^{1.75}-1.1231{\mathfrak{t}}^3+0.1070{\mathfrak{t}}^{3.25}-0.1876{\mathfrak{t}}^{3.5}+0.1667{\mathfrak{t}}^{3.75},\hfill \end{array}\right.$$

(67)

and

$$\left\{\begin{array}{cc}{\mu }_2\left(\mathfrak{t}\right)=\hfill & 0.1103{\mathfrak{t}}^{1.75}-0.1247{\mathfrak{t}}^3+0.1054{\mathfrak{t}}^{3.25}-0.0063{\mathfrak{t}}^{3.5}+0.0054{\mathfrak{t}}^{3.75}-0.0343{\mathfrak{t}}^4\hfill \\ \hfill & +0.0322{\mathfrak{t}}^{4.25}-0.0213{\mathfrak{t}}^{4.5}+0.0204{\mathfrak{t}}^{4.75}-0.0017{\mathfrak{t}}^5,\hfill \\ {\nu }_2\left(\mathfrak{t}\right)=\hfill & 0.7307{\mathfrak{t}}^{1.75}-1.1231{\mathfrak{t}}^3+0.1070{\mathfrak{t}}^{3.25}-0.1876{\mathfrak{t}}^{3.5}+0.1667{\mathfrak{t}}^{3.75}-0.0436{\mathfrak{t}}^4\hfill \\ \hfill & +0.0378{\mathfrak{t}}^{4.25}-0.0293{\mathfrak{t}}^{4.5}+0.0254{\mathfrak{t}}^{4.75}-0.0212{\mathfrak{t}}^5.\hfill \end{array}\right.$$

(68)

The results are illustrated graphically in Figure 1.

Note to the following example in which we change the coefficients of the boundary conditions.

Example 2.

We design the following system of coupled sequential Navier FBVPs with new boundary conditions.

$$\left\{\begin{array}{c}{\displaystyle {}^C{\mathfrak{D}}^{\frac{5}{4}}\left({}^C{\mathfrak{D}}^{\frac{7}{4}}\mu \right)\left(\mathfrak{t}\right)=\frac{e^{-\mathfrak{t}}\left(|\nu \left(\mathfrak{t}\right)|+|{}^C{\mathfrak{D}}^{\frac{7}{4}}\nu \left(\mathfrak{t}\right)|\right)}{\left(9+e^{-\mathfrak{t}}\right)\left(1+|\nu \left(\mathfrak{t}\right)|+|{}^C{\mathfrak{D}}^{\frac{7}{4}}\nu \left(\mathfrak{t}\right)|\right)},(\mathfrak{t}\in \mathcal{O}:=[0,1]),}\hfill \\ {\displaystyle {}^C{\mathfrak{D}}^{\frac{3}{2}}\left({}^C{\mathfrak{D}}^{\frac{5}{4}}\nu \right)\left(\mathfrak{t}\right)=\frac{|\mu \left(\mathfrak{t}\right)|+|{}^C{\mathfrak{D}}^{\frac{8}{5}}\mu \left(\mathfrak{t}\right)|}{\left(\pi +e^{-2\mathfrak{t}}\right)\left(1+|\mu \left(\mathfrak{t}\right)|+|{}^C{\mathfrak{D}}^{\frac{8}{5}}\mu \left(\mathfrak{t}\right)|\right)},(\mathfrak{t}\in \mathcal{O}:=[0,1]),}\hfill \\ \sqrt{2}\mu \left(0\right)=6\mu \left(1\right)=\sqrt{2}{}^C{\mathfrak{D}}^{\frac{7}{4}}\mu \left(0\right)=3{}^C{\mathfrak{D}}^{\frac{7}{4}}\mu \left(1\right)=0,\hfill \\ \pi \nu \left(0\right)=3\nu \left(1\right)=\pi {}^C{\mathfrak{D}}^{\frac{5}{4}}\nu \left(0\right)={}^C{\mathfrak{D}}^{\frac{5}{4}}\nu \left(1\right)=0.\hfill \end{array}\right.$$

(69)

Similarly to the previous example, from the recurrence relationship (65) truncated, respectively, at $\mathfrak{s}=15,$ $\mathfrak{s}=20,$ the approximate solutions $({\mu }_1\left(\mathfrak{t}\right),{\nu }_1\left(\mathfrak{t}\right)),$ $({\mu }_2\left(\mathfrak{t}\right)$, and ${\nu }_2\left(\mathfrak{t}\right))$ are obtained by

$$\left\{\begin{array}{c}{\mu }_1\left(\mathfrak{t}\right)=0.1012{\mathfrak{t}}^{1.75}-0.1333{\mathfrak{t}}^3+0.0527{\mathfrak{t}}^{3.25}-0.0173{\mathfrak{t}}^{3.5}+0.0027{\mathfrak{t}}^{3.75},\hfill \\ {\nu }_1\left(\mathfrak{t}\right)=0.5321{\mathfrak{t}}^{1.75}-0.1231{\mathfrak{t}}^3+0.0357{\mathfrak{t}}^{3.25}-0.3689{\mathfrak{t}}^{3.5}+0.0556{\mathfrak{t}}^{3.75},\hfill \end{array}\right.$$

(70)

and

$$\left\{\begin{array}{cc}{\mu }_2\left(\mathfrak{t}\right)=\hfill & 0.1012{\mathfrak{t}}^{1.75}-0.1333{\mathfrak{t}}^3+0.0527{\mathfrak{t}}^{3.25}-0.0173{\mathfrak{t}}^{3.5}+0.0027{\mathfrak{t}}^{3.75}-0.0124{\mathfrak{t}}^4\hfill \\ \hfill & +0.0161{\mathfrak{t}}^{4.25}-0.00284{\mathfrak{t}}^{4.5}+0.0102{\mathfrak{t}}^{4.75}-0.0083{\mathfrak{t}}^5,\hfill \\ {\nu }_2\left(\mathfrak{t}\right)=\hfill & 0.5321{\mathfrak{t}}^{1.75}-0.1231{\mathfrak{t}}^3+0.0357{\mathfrak{t}}^{3.25}-0.3689{\mathfrak{t}}^{3.5}+0.0556{\mathfrak{t}}^{3.75}-0.0139{\mathfrak{t}}^4\hfill \\ \hfill & +0.0226{\mathfrak{t}}^{4.25}-0.0898{\mathfrak{t}}^{4.5}+0.0085{\mathfrak{t}}^{4.75}-0.0106{\mathfrak{t}}^5.\hfill \end{array}\right.$$

(71)

The results are illustrated graphically in Figure 2.

7. Conclusions

In this paper, by means of a fixed point theorem due to Krasnoselskii, we studied the existence criterion for solutions of a system of coupled sequential Navier FBVPs and then investigated its uniqueness by terms of the contraction principle due to Banach. After that, the necessary criteria for H-U-type stability of solutions to such a system of coupled sequential FBVPs arising in Navier model (1) have been derived. Next, the approximate solutions of the given system were computed and founded via the generalized differential transformation, and we illustrated some results graphically in relation to the given system of Navier FBVPs in two different numerical examples based on it. If we compare the exact solution to the approximate solution in these examples, then we find that the generalized differential transform method yields the accurate outcomes for the assumed fractional nonlinear sequential Navier problem. The simplicity and accuracy of this numerical method reveal its applicability for approximating the solutions of different fractional systems arising in various real-world models. Even in future studies, we aim to conduct a similar research study for several nonlinear systems of real-world models in the framework of the fractional FBVPs supplemented with generalized operators having non-singular kernels .