Analytical Study of the Complexities in a Three Species Food Web Model with Modified Caputo–Fabrizio Operator

Abstract

This article presents the analytical study of the three species fractional food web model in the framework of the Modified Caputo–Fabrizio operator. With the help of fixed point theory, the existence and uniqueness results are investigated for the fractional order model. To obtain the approximate solution for the suggested model, the well-known Laplace–Adomian decomposition method is used. The solutions are validated through simulations with a variety of fractional orders and initial values, where the complex nature of the system can be observed. The technique used here can be easily used to study a range of complex problems in different branches of science. From the figures, it can be observed that, at integer higher fractional order, there are a number of oscillations in the system and the system behaves chaotically, while, at lower fractional orders, the oscillation amplitudes decrease, resulting in the faster converging towards the equilibrium point. According to the results, the Modified Caputo–Fabrizio fractional-order derivative may be used in a variety of future fractional dynamics scenarios.

1. Introduction

Fractional calculus (FC) is widely employed in a number of scientific fields. On the basis of this importance, many other definitions have been introduced in literature [1,2]. The non-local property of fractional derivatives (FD) has shown to be a key characteristic in a wide range of FC implementations [3,4,5]. The integer order derivative of a function at a given point can be approximated using nearby data; however, the FD mandates the entire history beginning at the origin. This non-local nature of the FD is significant in simulating the system’s storage and heredity features [6]. Thus, the models that use FD are more convincing than those that use integer order. The FD’s second advantage is its ability to imitate intermediary activities. The integer order derivatives cannot capture the realistic occurrences in different physical conditions, for instance the fluid motion in porous media, because these processes are intermediate. Due to the inability of analytical methodologies to solve much fractional-order nonlinear systems, a variety of methods have been introduced to obtain a reasonable solution to non-integer order phenomena [7,8,9,10,11,12].

In FC, a number of fractional order operators are available in the literature. From Caputo’s perspective, these operators are Caputo, Caputo–Fabrizio (CF), and Atangana–Baleanu (ABC) [13,14,15,16]. The CF operator was further improved by Caputo and Fabrizio [15], which is very valuable as compared to the CF operator presented earlier. The characteristic of this new fractional derivative operator was discussed in detail [17] These operators are very important since integer order operators are unable to analyse fractional order nonlinear differential equations (FNDEs) in order to obtain explicit solutions due to their complexity. Due to this disadvantage of classical operators, the coefficients of the series solutions of FNDEs must be obtained using the best possible numerical method [18,19,20].

The differential equations (DEs) have become an important branch of applied mathematics, which is used to explain a variety of physical phenomena. Complex systems are one of the most important and rapidly growing domains where DEs are used extensively. Every point in dimensional space can be characterised by the dynamical system over time. The status of the mathematical model at any given time t can be predicted if it is a real-world phenomena throughout in the state of a dynamical model. Fractional-order calculus has been found to be a very important for comprehending complicated dynamical models with nonlinear features in the study and the forecasting of such systems [21,22].

The predator–prey model is a significant class of ecological models. The predator–prey model is one of the most important issues in ecology because of the prevalent importance and persistence of the convoluted food webs networks. Research on dynamical behavior of prey–predator models became a major research subject for biologists and mathematicians since the works of Lotka and Volterra [23]. Various biological populations in our natural ecosystem employ strategies such as collective defense, refuging, fleeing, and so on to ensure their own existence or find nourishment. Many scientists pay close attention to the prey’s group defense mechanism in order to better understand the rule of predator–prey interaction. Falconi et al. [24] studied the stability and the global dynamics of a prey–predator model with group defense, and Raw et al. [25] investigated the complex dynamical behavior of a prey–predator model involving group defense.

Several features of the predator–prey models have been analysed and investigated by various researchers. For instance, Yahuz et al. studied the stability of the fractional order predator–prey system by considering the harvesting rate [26]. The predator–prey involving the social behavior of prey and infections in predator populations from the perspective of fractional operators are studied in [27]. Similarly, the chaotic dynamics, stability, bifurcations, and the influence of the prey escaping from the prey herd in fractional sense have been reported in the literature [28,29,30,31].

This study will demonstrate that some modifications to the definition of the CF fractional-order derivative are necessary to resolve some inconsistencies between the definition provided by Caputo and Fabrizio and the corresponding integral operator for the CF derivative. This study introduces a straightforward modification to the CF fractional-order derivative and demonstrates how this Modified-CF fractional-order derivative enables the use of various well-known techniques, including the Laplace–Adomian decomposition approach in order to obtain approximate analytical solutions to fractional differential equations when the new Modified-CF fractional-order derivative is taken into account. According to the simulations, the system exhibits a variety of oscillations and exhibits chaotic behavior at higher fractional orders, whereas, at lower fractional orders, the oscillations’ amplitudes diminish and the system rapidly approaches equilibrium.

2. Model Formulation

Mathematical formulations are used to describe the dynamics of a three-species food web model that includes cannibalism in the top predator in addition to a stage structure. ${\mathbb{W}}_1\left(t\right)$ shows the prey population, ${\mathbb{W}}_2\left(t\right)$ shows the density of the intermediate predator at time t and ${\mathbb{W}}_3\left(t\right)$, and ${\mathbb{W}}_4\left(t\right)$ stands for the mature and immature population of the species with high tropic level (top predator) at t.

Inspired from the above literature, a three species food web model is studied with a modified fractional-order CF operator. The Laplace–Adomian technique is used to obtain a semi analytical solution, and a flowchart for the model is presented in Figure 1. Here, consider the food web model [32]

$$\begin{array}{ccc}\hfill \dot{{\mathbb{W}}_1}& =& \gamma {\mathbb{W}}_1-\frac{\gamma {\mathbb{W}}_1^2}{\mathcal{L}}-{\rho }_1{\mathbb{W}}_1{\mathbb{W}}_2,\hfill \\ \hfill \dot{{\mathbb{W}}_2}& =& {\rho }_1{\mathrm{𝔟}}_1{\mathbb{W}}_1{\mathbb{W}}_2-{\rho }_2{\mathbb{W}}_2{\mathbb{W}}_3-{\omega }_1{\mathbb{W}}_2,\hfill \\ \hfill \dot{{\mathbb{W}}_3}& =& {\rho }_2{\mathrm{𝔟}}_2{\mathbb{W}}_2{\mathbb{W}}_3+{\rho }_3{\mathrm{𝔟}}_3{\mathbb{W}}_3{\mathbb{W}}_4+{\mathbb{W}}_4-{\omega }_2{\mathbb{W}}_3,\hfill \\ \hfill \dot{{\mathbb{W}}_4}& =& \alpha {\mathbb{W}}_3-{\mathbb{W}}_4-{\rho }_3{\mathbb{W}}_3{\mathbb{W}}_4-{\omega }_3{\mathbb{W}}_4,\hfill \end{array}$$

(1)

with the initial values ${\mathbb{W}}_1\left(0\right),{\mathbb{W}}_2\left(0\right),{\mathbb{W}}_3\left(0\right),{\mathbb{W}}_4\left(0\right)\ge 0$. In model (1), the used parameters are $\gamma$, which is the rate of birth, and carrying capacity $\mathcal{L}$; the prey population logistically grows. According to the functional response of Lotka–Volterra (LV), the ${\mathbb{W}}_2\left(t\right)$ consumes prey at a low level, with the maximum rate of attacks ${\rho }_1$ and the rate of conversion ${\mathrm{𝔟}}_1$. With no food source, it exponentially decays due to natural mortality rate ${\omega }_1$. Top predators are of two types, mature and immature. The immature population is considered to exponentially grow with their father and mother shown by the mature populace having the rate of growth $\alpha$, while some portion of them grows up to become adults with the rate of growth. Furthermore, the mature and immature populace face the natural rates of death with ${\omega }_2$ and ${\omega }_3$ accordingly. The maximum rate of attack is ${\rho }_2$ and the rate of conversion ${\mathrm{𝔟}}_2$; from the model, the mature top attacks on intermediate predator according to the functional response of LV. When there is a shortage of their preferred food, it cannibalises on the immature top predator based on the functional response LV with a maximum rate of attack ${\rho }_3$ and the rate of conversion ${\mathrm{𝔟}}_3$. In the fractional CF operator sense, the above model can be expressed as

$$\begin{array}{ccc}\hfill {}^{CF}{D}^{\mathrm{𝔡}}{\mathbb{W}}_1\left(t\right)& =& \gamma {\mathbb{W}}_1-\frac{\gamma {\mathbb{W}}_1^2}{\mathcal{L}}-{\rho }_1{\mathbb{W}}_1{\mathbb{W}}_2,\hfill \\ \hfill {}^{CF}{D}^{\mathrm{𝔡}}{\mathbb{W}}_2\left(t\right)& =& {\rho }_1{\mathrm{𝔟}}_1{\mathbb{W}}_1{\mathbb{W}}_2-{\rho }_2{\mathbb{W}}_2{\mathbb{W}}_3-{\omega }_1{\mathbb{W}}_2,\hfill \\ \hfill {}^{CF}{D}^{\mathrm{𝔡}}{\mathbb{W}}_3\left(t\right)& =& {\rho }_2{\mathrm{𝔟}}_2{\mathbb{W}}_2{\mathbb{W}}_3+{\rho }_3{\mathrm{𝔟}}_3{\mathbb{W}}_3{\mathbb{W}}_4+{\mathbb{W}}_4-{\omega }_2{\mathbb{W}}_3,\hfill \\ \hfill {}^{CF}{D}^{\mathrm{𝔡}}{\mathbb{W}}_4\left(t\right)& =& \alpha {\mathbb{W}}_3-{\mathbb{W}}_4-{\rho }_3{\mathbb{W}}_3{\mathbb{W}}_4-{\omega }_3{\mathbb{W}}_4,\hfill \end{array}$$

(2)

where $0<\varsigma \le 1$, and ${}^{CF}D$ represents the CF operator.

3. Equilibria Points and Stability Analysis

Here, equilibria points and stability analysis of the suggested model (1) are investigated. The following equilibrium points are observed for the system (1) as:

• The zero equilibrium point which always exists, $E_0=(0,0,0,0)$;
• The axial equilibrium points always exists, $E_1=(\mathcal{L},0,0,0)$;
• The top predator free equilibrium point $E_2=(\frac{{\omega }_1}{{\mathrm{𝔟}}_1{\rho }_1},\frac{\gamma }{{\rho }_1\mathcal{L}}(\mathcal{L}-\frac{{\omega }_1}{{\rho }_1{\mathrm{𝔟}}_1}),0,0)$ exists if $\mathcal{L}>\frac{{\omega }_1}{{\rho }_1{\mathrm{𝔟}}_1}$.${\mathbb{W}}_1,{\mathbb{W}}_2,{\mathbb{W}}_3,{\mathbb{W}}_4$
• The endemic equilibrium points are denoted by $E_3=({\mathbb{W}}_1^{\ast },{\mathbb{W}}_2^{\ast },{\mathbb{W}}_3^{\ast },{\mathbb{W}}_4^{\ast })$

  $$\begin{array}{ccc}\hfill {\mathbb{W}}_2^{\ast }& =& \frac{\gamma (\mathcal{L}-{\mathbb{W}}_1^{\ast })}{{\rho }_1\mathcal{L}}\hfill \end{array}$$

  (3)

  $$\begin{array}{ccc}\hfill {\mathbb{W}}_3^{\ast }& =& \frac{{\mathrm{𝔟}}_1{\rho }_1{\mathbb{W}}_1^{\ast }-{\omega }_1}{{\rho }_2}\hfill \end{array}$$

  (4)

  $$\begin{array}{ccc}\hfill {\mathbb{W}}_4^{\ast }& =& \frac{\phi ({\mathrm{𝔟}}_1{\rho }_1{\mathbb{W}}_1^{\ast }-{\omega }_1)}{{\rho }_2(+{\omega }_3)+{\rho }_3({\mathrm{𝔟}}_1{\rho }_1{\mathbb{W}}_1^{\ast }-{\omega }_1)}\hfill \end{array}$$

  (5)

where ${\mathbb{W}}_1^{\ast }$ is a positive root of the following third order polynomial equation:

$$\begin{array}{c}\hfill {\mathrm{D}}_1{\mathbb{W}}_1^{{\ast }^3}+{\mathrm{D}}_2{\mathbb{W}}_1^{{\ast }^2}+{\mathrm{D}}_3{\mathbb{W}}_1^{\ast }+{\mathrm{D}}_4=0,\end{array}$$

(6)

where ${\mathrm{D}}_1=-\gamma {\mathrm{𝔟}}_2{\mathrm{𝔟}}_1^2{\rho }_1^2{\rho }_2{\rho }_3<0$

$$\begin{array}{ccc}\hfill {\mathrm{D}}_2& =& {\mathrm{𝔟}}_1{\rho }_1[\gamma {\mathrm{𝔟}}_2{\rho }_2({\rho }_3{\mathrm{𝔟}}_3-{\rho }_2(+{\omega }_3))+{\mathrm{𝔟}}_1{\rho }_1^2{\rho }_3\mathcal{L}({\rho }_3\phi -{\omega }_2)+\gamma {\rho }_2{\rho }_3({\mathrm{𝔟}}_1{\rho }_1{\mathrm{𝔟}}_2\mathcal{L}+{\omega }_1)],\hfill \\ \hfill {\mathrm{D}}_3& =& \gamma {\rho }_2^2(+{\omega }_3)({\mathrm{𝔟}}_1{\rho }_1{\mathrm{𝔟}}_2\mathcal{L}+2{\omega }_1)-2\mathcal{L}{\mathrm{𝔟}}_1{\rho }_1{\rho }_3{\omega }_1(\gamma {\mathrm{𝔟}}_2{\rho }_2-{\omega }_2{\rho }_1)-\mathcal{L}{\mathrm{𝔟}}_1{\rho }_1^2{\rho }_2[{\omega }_2{\omega }_3+({\omega }_2-\phi )]\hfill \\ & -& {\mathrm{𝔟}}_1{\rho }_1^2{\omega }_1\phi \mathcal{L}({\mathrm{𝔟}}_3{\rho }_3+1)-\gamma {\rho }_2{\rho }_3{\omega }_1^2,\hfill \\ \hfill {\mathrm{D}}_4& =& \gamma \mathcal{L}{\mathrm{𝔟}}_2{\rho }_2{\mathrm{𝔟}}_1{({\rho }_3{\omega }_1-{\rho }_2(+{\omega }_3))}_{\mathcal{L}}{\rho }_1{\rho }_3{\omega }_1^2({\mathrm{𝔟}}_3\phi -{\omega }_2)+\mathcal{L}{\rho }_1{\omega }_1({\rho }_2{\omega }_2(+{\omega }_3)-\phi )\hfill \end{array}$$

Thus, from the above relations, its straightforward that $E_3$ exists if and only if the condition is satisfied:

$$\begin{array}{c}\hfill {\overline{\mathbb{W}}}_1<{\mathbb{W}}_1^{\ast }<\mathcal{L},\end{array}$$

(7)

with one from the below relation

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\hfill {\mathrm{D}}_2<0,\mathrm{and}{\mathrm{D}}_4>0,\\ \hfill {\mathrm{D}}_3>0,\mathrm{and}{\mathrm{D}}_4>0.\end{array}\right.\end{array}$$

(8)

For the local and global stability analysis, the author advised the reader to see [32] for more details.

4. Basic Results

Here, some basic definitions are presented from the literature [14,15,17].

Definition 1.

Consider a function $\mathrm{\Theta }\in {\mathrm{H}}^1[0,T]$, and the CF of fractional derivative for order $\mathrm{𝔡}\in (0,1),$ is presented as

$${}^{CF}{D}^{\mathrm{𝔡}}\left(\mathrm{\Theta }\left(t\right)\right)=\frac{\mathfrak{M}(\mathrm{𝔡})}{(1-\mathrm{𝔡})}{\int }_0^t{\mathrm{\Theta }}^{\prime }\left(t\right)exp\left[-\mathrm{𝔡}\frac{t-\rho }{1-\mathrm{𝔡}}\right]d\rho ,$$

where $\mathfrak{M}(\mathrm{𝔡})$ is $\mathfrak{M}(\mathrm{𝔡})=\frac{2}{2-\mathrm{𝔡}},0<\mathrm{𝔡}\le 1$. Moreover, $\mathfrak{M}\left(0\right)=1.$ If Θ∉ in ${\mathrm{H}}^1(0,T)$, then the operator of CF can be as

$${}^{CF}{D}^{\mathrm{𝔡}}(\mathrm{\Theta }\left(t\right))=\frac{\mathfrak{M}(\mathrm{𝔡})}{(1-\mathrm{𝔡})}{\int }_0^t(\mathrm{\Theta }\left(t\right)-\mathrm{\Theta }(\rho ))exp\left[-\mathrm{𝔡}\frac{t-\rho }{1-\mathrm{𝔡}}\right]d\rho .$$

Definition 2.

Consider a smooth function $F\left(t\right):[a,\infty )\to R$ with $a<0,t>0$, and the modified operator of CF can be defined as

$$\begin{array}{c}\hfill {}^{Modified-CF}{D}^{\mathrm{𝔡}}\left(F\left(t\right)\right)=\frac{1}{1-\mathrm{𝔡}}{\int }_0^t\left({\psi }_0^{\prime }\left(\xi \right)\beta (t-\xi )+F^{\prime }\left(\xi \right)exp\left(-\frac{\mathrm{𝔡}}{1-\mathrm{𝔡}}(t-\xi )\right)\right)d\xi \end{array}$$

(9)

where

$$\begin{array}{c}\hfill {\psi }_0^{\prime }\left(\xi \right)={\int }_a^0{F}^{\prime }\left(\xi \right)exp\left(-\frac{\mathrm{𝔡}}{1-\mathrm{𝔡}}(t-\xi )\right)d\xi \end{array}$$

(10)

Definition 3.

Consider $\mathrm{𝔡}\in ]0,1[;$ then, the integral of CF with order 𝔡 of Θ can be presented as

$${}^{CF}{I}^{\mathrm{𝔡}}[\mathrm{\Theta }\left(t\right)]=G\mathrm{\Theta }\left(t\right)+\overline{G}{\int }_0^t\mathrm{\Theta }(\rho )d\rho ,t\ge 0.$$

The classical integral for $\mathrm{\Theta }$ may be obtained when $\mathrm{𝔡}=1$; here,

$$G=\frac{(1-\mathrm{𝔡})}{\mathfrak{M}(\mathrm{𝔡})},\overline{G}=\frac{\mathrm{𝔡}}{\mathfrak{M}(\mathrm{𝔡})}.$$

Definition 4.

The transformation of Laplace for CF operator can be written as

$$\mathcal{L}\left\{{}^{Modified-CF}{D}_x^{\mathrm{𝔡}}F(x,t)\right\}=\frac{a(2-\mathrm{𝔡})}{2(s+(1-s)\mathrm{𝔡})}\left[{\xi }^{n+1}\overline{F}(s,\tau )-{\displaystyle \sum _{k=0}^n}{s}^{n-k}\left\{\frac{{\partial }^{k}F(0,t)}{\partial x^k}\right\}\right],$$

where $n=[\mathrm{𝔡}]+1$.

5. Theoretical Results

This part deals with the existence and uniqueness of solution for the considered model (2). Before proving the solution, some notations and lemma are presented.

Let us suppose $F=[0,1]$ and $\mathbb{C}\left(F\right)$ shows the space which contains continuous functions on $F$. Next, to suppose a Banach space, set $B=\mathsf{U}\left(t\right)/\mathsf{U}\left(t\right)\in \mathbb{C}\left(F\right)$ with ${\left|\right|\mathsf{U}\left(t\right)\left|\right|}_b\le {max}_{t\in F}\left|\mathsf{U}\left(t\right)\right|$. For clarity, one may suppose $\mathsf{U}\left(t\right)=\mathrm{\Psi }(t,\mathrm{\Omega }\left(t\right)),\nu \left(t\right)={\mathrm{\Psi }}_1(t,{\mathrm{\Omega }}_1\left(t\right))$ and $\mathsf{U}\left(0\right)=\mathrm{\Omega }\left(0\right)={\mathsf{U}}_0$ and $\mathcal{V}\left(0\right)={\mathrm{\Omega }}_1\left(0\right)={\mathcal{V}}_0$. The proposed system (2) in the form of integration is

$$\mathsf{U}\left(t\right)={\mathsf{U}}_0+\frac{2(1-\mathrm{𝔡})}{(2-\mathrm{𝔡})M(\mathrm{𝔡})}(\mathrm{\Psi }(t,\mathrm{\Omega }\left(t\right)))+\frac{(2\mathrm{𝔡})}{(2-\mathrm{𝔡})M(\mathrm{𝔡})}{\int }_0^t\mathrm{\Psi }(\wp ,\mathrm{\Omega }(\wp ))d\wp .$$

Here, to consider an operator, $\mathbf{T}:B\to B$ is defined as

$$\mathbf{T}\mathsf{U}\left(t\right)={\mathsf{U}}_0+\frac{2(1-\mathrm{𝔡})}{(2-\mathrm{𝔡})M(\mathrm{𝔡})}(\mathrm{\Psi }(t,\mathrm{\Omega }\left(t\right))+\frac{(2\mathrm{𝔡})}{(2-\mathrm{𝔡})M(\mathrm{𝔡})}{\int }_0^t\mathrm{\Psi }(\wp ,\mathrm{\Omega }(\wp ))d\wp ,$$

then, the operator $\mathbf{T}$ has the same fixed-point (FP) as (2).

Theorem 1.

Let us consider that a function $f:\mathrm{\Psi }\times ℛ\to ℛ$ must satisfy one from the below:

$\left({\mathscr{H}}_1\right)$

Consider a non-negative function $g\left(t\right)\in L[0,1]$ such that

$$|\mathrm{\Psi }(t,x)|\le h\left(t\right)+c_0{\left|x\right|}^{\wp },here,c_0\ge 0,0<\wp <1.$$

$\left({\mathscr{H}}_2\right)$

Consider a function Ψ satisfies $|\mathrm{\Psi }(t,x)|\le c_0{\left|x\right|}^{\wp }$, where $c_0>0$, $\wp >1$. Then, the proposed system (2) has at least one solution.

Proof. 

To obtain the required results, using the fixed point theorem of Schauder, first consider that the assumption $\left({\mathscr{H}}_1\right)$ is satisfied. Let us suppose $\mathcal{G}=\left\{\mathsf{U}\left(t\right)\right|\mathsf{U}\left(t\right)\in B,\parallel \mathsf{U}\left(t\right){\parallel }_B\le k,t\in F\}$, where $k\ge max{\left(2A{c}_0\right)}^{\left(\frac{1}{1-\wp }\right)},2l$ and $l={max}_{y\in F}\left({\mathsf{U}}_0+\frac{4(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)g\left(t\right)}+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t\left|g(\wp )\right|d\wp \right)$. Clearly, in B, $𝓖$ is a ball. Next, show that $\mathbf{T}:𝓖\to 𝓖$.

∀$u\in G$ to obtain

$$\begin{array}{cc}& \left|\mathbf{T}\mathsf{U}\left(t\right)\right|=|{\mathsf{U}}_0+\frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}(\mathrm{\Psi }(t,\mathrm{\Omega }\left(t\right)))+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t\mathrm{\Psi }(\wp ,\mathrm{\Omega }(\wp ))d\wp |\hfill \\ & \le {\mathsf{U}}_0+\frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}|\mathrm{\Psi }(t,\mathrm{\Omega }\left(t\right))|+\frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t\mathrm{\Psi }(\wp ,\mathrm{\Omega }(\wp ))d\wp \hfill \\ & \le {\mathsf{U}}_0+\frac{4(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}\left(g\left(t\right)+c_0{k}^{\wp }\right)+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t(g(\wp )+c_0{k}^{\wp })d\wp \hfill \\ & \le {\mathsf{U}}_0+\frac{4(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}\left(g\left(t\right)+c_0{k}^{\wp }\right)+\frac{2\mathbf{𝔡}{c}_0{k}^{\wp }t}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t(g(\wp )d\wp \hfill \\ & \le {\mathsf{U}}_0+\frac{4(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}g\left(t\right)+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t(g(\wp )d\wp +\left(\frac{4(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}+\frac{2\mathbf{𝔡}t}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}\right)c_0{k}^{\wp }\hfill \\ & \le {\mathsf{U}}_0+\frac{4(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}g\left(t\right)+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^{t}g(\wp )d\wp +\left(\frac{4(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}+\frac{2c_0{k}^{\wp }}{M\left(\mathbf{𝔡}\right)}\right).\hfill \end{array}$$

Therefore,

$$\begin{array}{cc}\hfill & {\left|\right|\mathbf{T}\mathsf{U}\left(t\right)\left|\right|}_B=\underset{t\in F}{max}\left|\mathbf{T}\mathsf{U}\left(t\right)\right|\hfill \\ & \le l+\frac{2c_0{k}^{\wp }}{M\left(\mathbf{𝔡}\right)}=l+A{c}_0{k}^{\wp }\le \frac{k}{2}+\frac{k}{2}=k.\hfill \end{array}$$

Thus, the operator $\mathbf{T}\mathsf{U}\left(t\right)$ is continuous on $F$.

Now, to consider that the assumption $\left({\mathscr{H}}_2\right)$ is also satisfied, choose $0\le k\le \left(\frac{1}{A{c}_0}\right)(\frac{1}{\wp -1})$. In the same fashion, by repeating the above procedure, follow:

$${\left|\right|\mathbf{T}\mathsf{U}\left(t\right)\left|\right|}_B\le A{c}_0{k}^{\wp }\le k.$$

As a result, to obtain $\mathbf{T}:G\to G$, clearly the operator $\mathbf{T}$ is continuous because of the continuity of $\mathrm{\Psi }$.

Next, show that the operator T is completely continuous. Let $R={max}_{t\in F}\left|F(t,\mathrm{\Omega }\left(t\right))\right|$, for any $F\in G$. Let $t_1,t_2\in F$ such that $t_1<t_2$.

In addition, let ${\mathrm{\Psi }}_1=\frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}$ and ${\mathrm{\Psi }}_2=\frac{2\mathbf{𝔡}}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}$,

$$\begin{array}{cc}\hfill & \mid T\mathrm{\Psi }\left(t_2\right)-T\mathrm{\Psi }\left(t_1\right)\mid =\mid {\mathsf{U}}_0+{\mathrm{\Psi }}_1\left[\mathsf{U}(t_2,\mathrm{\Psi }\left(t_2\right))\right]+{\mathrm{\Psi }}_2{\int }_0^{t_2}f(\wp ,u(\wp ))d\wp \hfill \\ & -{\mathsf{U}}_0-{\mathrm{\Psi }}_1\left[\mathsf{U}(t_1,\mathrm{\Psi }\left(t_1\right))\right]+{\mathrm{\Psi }}_2{\int }_0^{t_2}\mathsf{U}(\wp ,u(\wp ))d\wp \mid \hfill \\ & =\mid {\mathrm{\Psi }}_1[\mathsf{U}(t_2,\mathrm{\Psi }\left(t_2\right))-\mathsf{U}(t_1,\mathrm{\Psi }\left(t_1\right))]+{\mathrm{\Psi }}_2{\int }_{t_1}^{t_2}\mathrm{\Psi }(\wp ,u(\wp ))d\wp \mid \hfill \\ & \le {\mathrm{\Psi }}_1\mid \mathsf{U}(t_2,\mathrm{\Psi }\left(t_2\right))\mid +{\mathrm{\Psi }}_1\mid \mathsf{U}(t_1,\mathrm{\Psi }\left(t_1\right))\mid +{\mathrm{\Psi }}_2{\int }_{t_1}^{t_2}\mid \mathsf{U}(\wp ,u(\wp ))\mid d\wp \hfill \\ & \le 2R{\mathrm{\Psi }}_1+R{\mathrm{\Psi }}_2{\int }_{t_1}^{t_2}d\wp =R(2{\mathrm{\Psi }}_1+{\mathrm{\Psi }}_2(t_2-t_1)).\hfill \end{array}$$

Use the uniform continuity of a function $(t_2-t_1)$ on interval $F$, to obtain that the operator $\mathbf{T}\mathcal{G}$ is equi-continuous. From this, the function is also uniformly bounded as $\mathbf{T}𝒢\subseteq \mathcal{G}$; therefore, the operator $\mathbf{T}$ is completely continuous. Hence, by utilizing the fixed point theorem of Schauder, a solution of Equation (2) ∃ in the set $\mathcal{G}$. □

Corollary 1.

Supposing that $\mathsf{U}$ is a continuous bounded function on $F\times \mathrm{R}$, then the proposed problem (2) has at least one solution.

Proof. 

We know that the function $\mathsf{U}$ is continuous and bounded on $F\times \mathrm{R}$, and there exists $L>0$, which satisfies the inequality $\left|\mathsf{U}\right|<L$. Suppose $h\left(t\right)=L,c_0=0$ in $\left({\mathscr{H}}_1\right)$ of Theorem 1; then, the proposed problem (2) has at least one solution.

In the next theorem, apply the Banach contraction approach to prove the uniqueness of solution for (2). □

Theorem 2.

Consider that a function $\mathsf{U}:F\times \mathrm{R}\to \mathrm{R}$ is continuous, and satisfies the condition below:

$\left({\mathscr{H}}_3\right)$

Suppose a function $h\left(t\right)\in L[0,1]$ exists which is positive, ∋

$$\mid \mathsf{U}\left(t\right)\mid \le h\mid \mathrm{\Psi }(t,\mathrm{\Omega }(t\left)\right)\mid ,0\le t\le 1,$$

and also a function $\mathsf{U}$ satisfying $\mathsf{U}\left(0\right)=0$

$\left(H_{\mathscr{H}}\right)$

Suppose that $\wp ={max}_{t\in F}\mid \frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}h\left(t\right)+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t\left|h(\wp )\right|d\wp \mid <1$; then, the proposed problem (2) solution is unique.

Proof. 

Before starting the proof, define the operator $\mathbf{T}$ as

$$\begin{array}{cc}\hfill & \mathbf{T}\mathsf{U}\left(t\right)={\mathsf{U}}_0+\frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}\mathrm{\Psi }(t,\mathrm{\Omega }\left(t\right))+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t|\mathrm{\Psi }(\wp ,u(\wp ))|d\wp ,\hfill \end{array}$$

For $\mathsf{U}\left(t\right)\in B$, as

$$\begin{array}{cc}\hfill & \mid \mathbf{T}\mathsf{U}\left(t\right)\mid =\mid {\mathsf{U}}_0+\frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}\mathrm{\Psi }(t,\mathrm{\Omega }\left(t\right))+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t\mid \mathrm{\Psi }(\wp ,u(\wp ))\mid d\wp ,\hfill \\ & \le \mid F\mid +\mid \frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}\mathrm{\Psi }(t,\mathrm{\Omega }\left(t\right))\mid +\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t\mid \mathrm{\Psi }(\wp ,u(\wp ))\mid d\wp \hfill \\ & \le \mid {\mathsf{U}}_0\mid +\frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}h\left(t\right)\mid \mathsf{U}\left(t\right)\mid +\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^{t}h(\wp )\mid u(\wp )d\wp \hfill \\ & \le \mid {\mathsf{U}}_0\mid +\left(\frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}h\left(t\right)+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^{t}h(\wp )d\wp \right)\parallel u\parallel ,\hfill \end{array}$$

one may have

$$\begin{array}{cc}\hfill & {\Vert \mathbf{T}\mathsf{U}\left(t\right)\Vert }_B\le \mid {\mathsf{U}}_0\mid +\left(\frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}h\left(t\right)+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^{t}h(\wp )d\wp \right)\parallel u\parallel \hfill \\ & \le \mid F\mid +\wp \Vert u\Vert \le \Vert u\Vert .\hfill \end{array}$$

Let $\mathsf{U}\left(t\right),v\left(t\right)\in B$; then,

$$\begin{array}{cc}\hfill & \mid \mathbf{T}\mathsf{U}\left(t\right)-\mathbf{T}\nu \left(t\right)\mid =\mid {\mathsf{U}}_0+\frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}\mathrm{\Psi }(t,\mathrm{\Omega }\left(t\right))+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t|\mathrm{\Psi }(\wp ,\mathrm{\Omega }(\wp ))|d\wp \hfill \\ & -{\mathcal{V}}_0-\frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\mathrm{\Psi }}_1(t,{\mathrm{\Omega }}_1\left(t\right))-\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t|{\mathrm{\Psi }}_1(\wp ,{\mathrm{\Omega }}_1(\wp ))d\wp \mid \hfill \\ & \le \frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}\mid \mathrm{\Psi }(t,\mathrm{\Omega }\left(t\right))-{\mathrm{\Psi }}_1(\wp ,{\mathrm{\Omega }}_1(\wp ))\mid +\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t\mid \mathrm{\Psi }(\wp ,\mathrm{\Omega }(\wp ))-{\mathrm{\Psi }}_1(\wp ,{\mathrm{\Omega }}_1(\wp ))d\wp \mid \hfill \\ & \le \frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}\mid \mathsf{U}\left(t\right)-\nu \left(t\right)\mid +\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t\mid \mathrm{\Omega }(\wp )-{\mathrm{\Omega }}_1(\wp )\mid d\wp \hfill \\ & \le \left(\frac{2(1-\mathbf{𝔡})}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}h\left(t\right)+\frac{\left(2\mathbf{𝔡}\right)}{(2-\mathbf{𝔡})M\left(\mathbf{𝔡}\right)}{\int }_0^t\mid h(\wp )d\wp \right)\mid \mathrm{\Omega }(\wp )-{\mathrm{\Omega }}_1(\wp )\mid \hfill \\ & \le \mathrm{\Psi }\parallel \mathrm{\Omega }(\wp )-\nu (\wp )\parallel \le \parallel \mathrm{\Omega }(\wp )-{\mathrm{\Omega }}_1(\wp )\parallel \hfill \end{array}$$

In view of $\wp <1$, from the above results, one may observe that the operator $\mathbf{T}$ is a contraction. Finally, from the Banach contraction principle, it shows that the operator $\mathbf{T}$ has only one fixed point. □

6. Numerical Solution

In this part of the paper, use the Laplace transform on both sides of the proposed model (2) to investigate the analytical results, as

$$\begin{array}{c}\hfill \begin{array}{cc}& \mathcal{L}{(}^{CF}{\mathrm{D}}_t^{\mathbf{𝔡}}{\mathbb{W}}_1\left(t\right))=\mathcal{L}\left(\gamma {\mathbb{W}}_1-\frac{\gamma {\mathbb{W}}_1^2}{\mathcal{L}}-{\rho }_1{\mathbb{W}}_1{\mathbb{W}}_2\right),\hfill \\ & \mathcal{L}{(}^{CF}{\mathrm{D}}_t^{\mathbf{𝔡}}{\mathbb{W}}_2\left(t\right))=\mathcal{L}\left({\rho }_1{\mathrm{𝔟}}_1{\mathbb{W}}_1{\mathbb{W}}_2-{\rho }_2{\mathbb{W}}_2{\mathbb{W}}_3-{\omega }_1{\mathbb{W}}_2\right),\hfill \\ & \mathcal{L}{(}^{CF}{\mathrm{D}}_t^{\mathbf{𝔡}}{\mathbb{W}}_3\left(t\right))=\mathcal{L}\left({\rho }_2{\mathrm{𝔟}}_2{\mathbb{W}}_2{\mathbb{W}}_3+{\rho }_3{\mathrm{𝔟}}_3{\mathbb{W}}_3{\mathbb{W}}_4+\wp {\mathbb{W}}_4-{\omega }_2{\mathbb{W}}_3\right),\hfill \\ & \mathcal{L}{(}^{CF}{\mathrm{D}}_t^{\mathbf{𝔡}}{\mathbb{W}}_4\left(t\right))=\mathcal{L}\left(\alpha {\mathbb{W}}_3-\wp {\mathbb{W}}_4-{\rho }_3{\mathbb{W}}_3{\mathbb{W}}_4-{\omega }_3{\mathbb{W}}_4\right).\hfill \end{array}\end{array}$$

(11)

Using the initial values, (2) yields

$$\begin{array}{c}\hfill \begin{array}{cc}& \mathcal{L}{(}^{CF}{\mathrm{D}}_t^{\mathbf{𝔡}}{\mathbb{W}}_1\left(t\right))=\frac{{\mathbb{W}}_1\left(0\right)}{s}+\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left(\gamma {\mathbb{W}}_1-\frac{\gamma {\mathbb{W}}_1^2}{\mathcal{L}}-{\rho }_1{\mathbb{W}}_1{\mathbb{W}}_2\right),\hfill \\ & \mathcal{L}{(}^{CF}{\mathrm{D}}_t^{\mathbf{𝔡}}{\mathbb{W}}_2\left(t\right))=\frac{{\mathbb{W}}_2\left(0\right)}{s}+\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left({\rho }_1{\mathrm{𝔟}}_1{\mathbb{W}}_1{\mathbb{W}}_2-{\rho }_2{\mathbb{W}}_2{\mathbb{W}}_3-{\omega }_1{\mathbb{W}}_2\right),\hfill \\ & \mathcal{L}{(}^{CF}{\mathrm{D}}_t^{\mathbf{𝔡}}{\mathbb{W}}_3\left(t\right))=\frac{{\mathbb{W}}_3\left(0\right)}{s}+\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left({\rho }_2{\mathrm{𝔟}}_2{\mathbb{W}}_2{\mathbb{W}}_3+{\rho }_3{\mathrm{𝔟}}_3{\mathbb{W}}_3{\mathbb{W}}_4+\wp {\mathbb{W}}_4-{\omega }_2{\mathbb{W}}_3\right),\hfill \\ & \mathcal{L}{(}^{CF}{\mathrm{D}}_t^{\mathbf{𝔡}}{\mathbb{W}}_4\left(t\right))=\frac{{\mathbb{W}}_4\left(0\right)}{s}+\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left(\alpha {\mathbb{W}}_3-\wp {\mathbb{W}}_4-{\rho }_3{\mathbb{W}}_3{\mathbb{W}}_4-{\omega }_3{\mathbb{W}}_4\right).\hfill \end{array}\end{array}$$

(12)

For the finite series form, consider the following:

$$\begin{array}{c}\hfill {\mathbb{W}}_1\left(t\right)={\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{1n}\left(t\right),{\mathbb{W}}_2\left(t\right)={\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{2n}\left(t\right),\\ \hfill {\mathbb{W}}_3\left(t\right)={\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{1n}\left(t\right),{\mathbb{W}}_4\left(t\right)={\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{4n}\left(t\right).\end{array}$$

(13)

For the nonlinear terms ${\mathbb{W}}_1^2$, apply the technique of Adomian polynomials as

$$\begin{array}{c}\hfill {\mathrm{A}}_{\upsilon }({\mathbb{W}}_1,{\mathbb{W}}_2)=\frac{1}{\upsilon !}\frac{d^{\upsilon }}{d{\lambda }^{\upsilon }}{\left[{\displaystyle \sum _{\kappa =0}^{\upsilon }}{\left({\lambda }^{\kappa }{\mathbb{W}}_{1\kappa }\right)}^2\right]}_{\lambda =0}.\end{array}$$

(14)

Several terms are computed by the Adomian polynomials, which are given as

$$\begin{array}{c}\hfill \begin{array}{cc}& \upsilon =0:{\mathrm{A}}_0(\left({\mathbb{W}}_1\right)={\mathbb{W}}_{10}^2,\hfill \\ & \upsilon =1:{\mathrm{A}}_1(\left({\mathbb{W}}_1\right)=2{\mathbb{W}}_{10}{\mathbb{W}}_{11},\hfill \\ & \upsilon =2:{\mathrm{A}}_2(\left({\mathbb{W}}_1\right)={\mathbb{W}}_{10}^2+2{\mathbb{W}}_{10}{\mathbb{W}}_{12},\hfill \\ & \upsilon =3:{\mathrm{A}}_3(\left({\mathbb{W}}_1\right)=2{\mathbb{W}}_{10}{\mathbb{W}}_{13}+2{\mathbb{W}}_{11}{\mathbb{W}}_{12},\hfill \\ & \upsilon =4:{\mathrm{A}}_4(\left({\mathbb{W}}_1\right)=2{\mathbb{W}}_{10}{\mathbb{W}}_{14}+2{\mathbb{W}}_{11}{\mathbb{W}}_{13}+{\mathbb{W}}_{12}^2,\hfill \end{array}\end{array}$$

(15)

and so on.

In addition, to express the nonlinear terms ${\mathbb{W}}_1{\mathbb{W}}_2$, one will also apply the Adomian polynomials as

$$\begin{array}{c}\hfill {\mathrm{B}}_{\upsilon }({\mathbb{W}}_1,{\mathbb{W}}_2)=\frac{1}{\upsilon !}\frac{d^{\upsilon }}{d{\lambda }^{\upsilon }}{\left[{\displaystyle \sum _{\kappa =0}^{\upsilon }}\left({\lambda }^{\kappa }{\mathbb{W}}_{1\kappa }\right)\left({\lambda }^{\kappa }{\mathbb{W}}_{2\kappa }\right)\right]}_{\lambda =0}.\end{array}$$

(16)

Several computed terms are below, which are obtained by using the Adomian polynomials technique:

$$\begin{array}{ccc}\hfill & & \upsilon =0:{\mathrm{B}}_0\left(({\mathbb{W}}_1,{\mathbb{W}}_2)\right)={\mathbb{W}}_{10}{\mathbb{W}}_{20},\hfill \\ & & \upsilon =1:{\mathrm{B}}_1\left(({\mathbb{W}}_1,{\mathbb{W}}_2)\right)={\mathbb{W}}_{10}{\mathbb{W}}_{21}+{\mathbb{W}}_{11}{\mathbb{W}}_{20},\hfill \\ & & \upsilon =2:{\mathrm{B}}_2\left(({\mathbb{W}}_1,{\mathbb{W}}_2)\right)={\mathbb{W}}_{10}{\mathbb{W}}_{22}+{\mathbb{W}}_{11}{\mathbb{W}}_{21}+{\mathbb{W}}_{12}{\mathbb{W}}_{20},\hfill \\ & & \upsilon =3:{\mathrm{B}}_3\left(({\mathbb{W}}_1,{\mathbb{W}}_2)\right)={\mathbb{W}}_{10}{\mathbb{W}}_{23}+{\mathbb{W}}_{11}{\mathbb{W}}_{22}+{\mathbb{W}}_{12}{\mathbb{W}}_{21}+{\mathbb{W}}_{13}{\mathbb{W}}_{20},\hfill \\ & & \upsilon =4:{\mathrm{B}}_4\left(({\mathbb{W}}_1,{\mathbb{W}}_2)\right)={\mathbb{W}}_{10}{\mathbb{W}}_{24}+{\mathbb{W}}_{11}{\mathbb{W}}_{23}+{\mathbb{W}}_{12}{\mathbb{W}}_{22}+{\mathbb{W}}_{13}{\mathbb{W}}_{21}+{\mathbb{W}}_{14}{\mathbb{W}}_{20},\hfill \end{array}$$

(17)

and so on. Similarly, for the nonlinear term ${\mathbb{W}}_2{\mathbb{W}}_3$, consider

$$\begin{array}{c}\hfill {\mathrm{C}}_{\upsilon }({\mathbb{W}}_2,{\mathbb{W}}_3)=\frac{1}{\upsilon !}\frac{d^{\upsilon }}{d{\lambda }^{\upsilon }}{\left[{\displaystyle \sum _{\kappa =0}^{\upsilon }}\left({\lambda }^{\kappa }{\mathbb{W}}_{2\kappa }\right)\left({\lambda }^{\kappa }{\mathbb{W}}_{3\kappa }\right)\right]}_{\lambda =0}.\end{array}$$

(18)

Several computed terms are below, which are obtained by using the Adomian polynomials technique:

$$\begin{array}{ccc}\hfill & & \upsilon =0:{\mathrm{C}}_0\left(({\mathbb{W}}_2,{\mathbb{W}}_3)\right)={\mathbb{W}}_{20}{\mathbb{W}}_{30},\hfill \\ & & \upsilon =1:{\mathrm{C}}_1\left(({\mathbb{W}}_2,{\mathbb{W}}_3)\right)={\mathbb{W}}_{20}{\mathbb{W}}_{31}+{\mathbb{W}}_{21}{\mathbb{W}}_{30},\hfill \\ & & \upsilon =2:{\mathrm{C}}_2\left(({\mathbb{W}}_2,{\mathbb{W}}_3)\right)={\mathbb{W}}_{20}{\mathbb{W}}_{32}+{\mathbb{W}}_{21}{\mathbb{W}}_{31}+{\mathbb{W}}_{22}{\mathbb{W}}_{30},\hfill \\ & & \upsilon =3:{\mathrm{C}}_3\left(({\mathbb{W}}_2,{\mathbb{W}}_3)\right)={\mathbb{W}}_{20}{\mathbb{W}}_{33}+{\mathbb{W}}_{21}{\mathbb{W}}_{32}+{\mathbb{W}}_{22}{\mathbb{W}}_{31}+{\mathbb{W}}_{23}{\mathbb{W}}_{30},\hfill \\ & & \upsilon =4:{\mathrm{C}}_4\left(({\mathbb{W}}_2,{\mathbb{W}}_3)\right)={\mathbb{W}}_{20}{\mathbb{W}}_{34}+{\mathbb{W}}_{21}{\mathbb{W}}_{33}+{\mathbb{W}}_{22}{\mathbb{W}}_{32}+{\mathbb{W}}_{23}{\mathbb{W}}_{31}+{\mathbb{W}}_{24}{\mathbb{W}}_{30},\hfill \end{array}$$

(19)

and so on.

In addition, for the nonlinear term ${\mathbb{W}}_3{\mathbb{W}}_4$, use the Adomian polynomials as

$$\begin{array}{c}\hfill {\mathrm{D}}_{\upsilon }({\mathbb{W}}_3,{\mathbb{W}}_4)=\frac{1}{\upsilon !}\frac{d^{\upsilon }}{d{\lambda }^{\upsilon }}{\left[{\displaystyle \sum _{\kappa =0}^{\upsilon }}\left({\lambda }^{\kappa }{\mathbb{W}}_{3\kappa }\right)\left({\lambda }^{\kappa }{\mathbb{W}}_{4\kappa }\right)\right]}_{\lambda =0}.\end{array}$$

(20)

Several computed terms are below, which are obtained by using the Adomian polynomials technique

$$\begin{array}{c}\hfill \begin{array}{cc}& \upsilon =0:{\mathrm{D}}_0\left(({\mathbb{W}}_3,{\mathbb{W}}_4)\right)={\mathbb{W}}_{30}{\mathbb{W}}_{40},\hfill \\ & \upsilon =1:{\mathrm{D}}_1\left(({\mathbb{W}}_3,{\mathbb{W}}_4)\right)={\mathbb{W}}_{30}{\mathbb{W}}_{41}+{\mathbb{W}}_{31}{\mathbb{W}}_{40},\hfill \\ & \upsilon =2:{\mathrm{D}}_2\left(({\mathbb{W}}_3,{\mathbb{W}}_4)\right)={\mathbb{W}}_{30}{\mathbb{W}}_{42}+{\mathbb{W}}_{31}{\mathbb{W}}_{41}+{\mathbb{W}}_{32}{\mathbb{W}}_{40},\hfill \\ & \upsilon =3:{\mathrm{D}}_3\left(({\mathbb{W}}_3,{\mathbb{W}}_4)\right)={\mathbb{W}}_{30}{\mathbb{W}}_{43}+{\mathbb{W}}_{31}{\mathbb{W}}_{42}+{\mathbb{W}}_{32}{\mathbb{W}}_{41}+{\mathbb{W}}_{33}{\mathbb{W}}_{40},\hfill \\ & \upsilon =4:{\mathrm{D}}_4\left(({\mathbb{W}}_3,{\mathbb{W}}_4)\right)={\mathbb{W}}_{30}{\mathbb{W}}_{44}+{\mathbb{W}}_{31}{\mathbb{W}}_{43}+{\mathbb{W}}_{32}{\mathbb{W}}_{42}+{\mathbb{W}}_{33}{\mathbb{W}}_{41}+{\mathbb{W}}_{34}{\mathbb{W}}_{40},\hfill \end{array}\end{array}$$

(21)

and so on. Now, putting (13), (15), (17), (19) and (21) in (12),

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{L}\left({}^{CF}{\mathrm{D}}_t^{\mathbf{𝔡}}{\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{1\upsilon }\left(t\right)\right)& =\frac{{\mathbb{W}}_1\left(0\right)}{s}+\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left(\gamma {\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{1\upsilon }\left(t\right)\right.\hfill \\ & \left.-\frac{\gamma {\left({\sum }_{\upsilon =0}^{\infty }{\mathrm{A}}_{\upsilon }\left({\mathbb{W}}_1\right)\right)}^2}{\mathcal{L}}-{\rho }_1\left({\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathrm{B}}_{\upsilon }({\mathbb{W}}_1,{\mathbb{W}}_2)\right)\right),\hfill \\ \hfill \mathcal{L}\left({}^{CF}{\mathrm{D}}_t^{\mathbf{𝔡}}{\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{2\upsilon }\left(t\right)\right)& =\frac{{\mathbb{W}}_2\left(0\right)}{s}+\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left({\rho }_1{\mathrm{𝔟}}_1\left({\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathrm{B}}_{\upsilon }({\mathbb{W}}_1,{\mathbb{W}}_2)\right)\right.\hfill \\ & \left.-{\rho }_2\left({\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathrm{C}}_{\upsilon }({\mathbb{W}}_2,{\mathbb{W}}_3)\right)-{\omega }_1{\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{2\upsilon }\right),\hfill \\ \hfill \mathcal{L}\left({}^{CF}{\mathrm{D}}_t^{\mathbf{𝔡}}{\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{3\upsilon }\left(t\right)\right)& =\frac{{\mathbb{W}}_3\left(0\right)}{s}+\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left({\rho }_2{\mathrm{𝔟}}_2\left({\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathrm{C}}_{\upsilon }({\mathbb{W}}_2,{\mathbb{W}}_3)\right)\right.\hfill \\ & \left.+{\rho }_3{\mathrm{𝔟}}_3\left({\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathrm{D}}_{\upsilon }({\mathbb{W}}_3,{\mathbb{W}}_4)\right)+\wp {\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{4\upsilon }-{\omega }_2{\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{3\upsilon }\right),\hfill \\ \hfill \mathcal{L}\left({}^{CF}{\mathrm{D}}_t^{\mathbf{𝔡}}{\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{4\upsilon }\left(t\right)\right)& =\frac{{\mathbb{W}}_4\left(0\right)}{s}+\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left(\wp {\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{3\upsilon }-\wp {\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{4\upsilon }\right.\hfill \\ & \left.-{\rho }_3\left({\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathrm{D}}_{\upsilon }({\mathbb{W}}_3,{\mathbb{W}}_4)\right)-{\omega }_3{\displaystyle \sum _{\upsilon =0}^{\infty }}{\mathbb{W}}_{4\upsilon }\right).\hfill \end{array}\end{array}$$

Now, compare similar terms from both sides as:

$$\begin{array}{ccc}\hfill & & \mathcal{L}{\mathbb{W}}_{10}\left(t\right)=\frac{{\mathbb{W}}_1\left(0\right)}{s},\mathcal{L}{\mathbb{W}}_{20}\left(t\right)=\frac{{\mathbb{W}}_2\left(0\right)}{s},\mathcal{L}{\mathbb{W}}_{30}\left(t\right)=\frac{{\mathbb{W}}_3\left(0\right)}{s},\mathcal{L}{\mathbb{W}}_{40}\left(t\right)=\frac{{\mathbb{W}}_4\left(0\right)}{s},\hfill \\ & & \mathcal{L}{\mathbb{W}}_{11}\left(t\right)=\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left(\gamma {\mathbb{W}}_{10\upsilon }-\frac{\gamma {\left({\mathrm{A}}_0\left({\mathbb{W}}_1\right)\right)}^2}{\mathcal{L}}-{\rho }_1{\mathrm{B}}_0({\mathbb{W}}_1,{\mathbb{W}}_2)\right),\hfill \\ & & \mathcal{L}{\mathbb{W}}_{21}\left(t\right)=\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left({\rho }_1{\mathrm{𝔟}}_1{\mathrm{B}}_0({\mathbb{W}}_1,{\mathbb{W}}_2)-{\rho }_2{\mathrm{C}}_0({\mathbb{W}}_2,{\mathbb{W}}_3)-{\omega }_1{\mathbb{W}}_{20}\right),\hfill \\ & & \mathcal{L}{\mathbb{W}}_{31}\left(t\right)=\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left({\rho }_2{\mathrm{𝔟}}_2{\mathrm{C}}_0({\mathbb{W}}_2,{\mathbb{W}}_3)+{\rho }_3{\mathrm{𝔟}}_3{\mathrm{D}}_0({\mathbb{W}}_3,{\mathbb{W}}_4)+\wp {\mathbb{W}}_{40}-{\omega }_2{\mathbb{W}}_{30}\right),\hfill \\ & & \mathcal{L}{\mathbb{W}}_{41}\left(t\right)=\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left(\alpha {\mathbb{W}}_{30}-\wp {\mathbb{W}}_{40}-{\rho }_3{\mathrm{D}}_0({\mathbb{W}}_3,{\mathbb{W}}_4)-{\omega }_3{\mathbb{W}}_{40}\right),\hfill \\ & & \mathcal{L}{\mathbb{W}}_{12}\left(t\right)=\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left(\gamma {\mathbb{W}}_{11}-\frac{\gamma {\left({\mathrm{A}}_1\left({\mathbb{W}}_1\right)\right)}^2}{\mathcal{L}}-{\rho }_1{\mathrm{B}}_1({\mathbb{W}}_1,{\mathbb{W}}_2)\right),\hfill \\ & & \mathcal{L}{\mathbb{W}}_{22}\left(t\right)=\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left({\rho }_1{\mathrm{𝔟}}_1{\mathrm{B}}_1({\mathbb{W}}_1,{\mathbb{W}}_2)-{\rho }_2{\mathrm{C}}_1({\mathbb{W}}_2,{\mathbb{W}}_3)-{\omega }_1{\mathbb{W}}_{21}\right),\hfill \\ & & \mathcal{L}{\mathbb{W}}_{32}\left(t\right)=\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left({\rho }_2{\mathrm{𝔟}}_2{\mathrm{C}}_1({\mathbb{W}}_2,{\mathbb{W}}_3)+{\rho }_3{\mathrm{𝔟}}_3{\mathrm{D}}_1({\mathbb{W}}_3,{\mathbb{W}}_4)+\wp {\mathbb{W}}_{41}-{\omega }_2{\mathbb{W}}_{31}\right),\hfill \\ & & \mathcal{L}{\mathbb{W}}_{42}\left(t\right)=\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left(\alpha {\mathbb{W}}_{31}-\wp {\mathbb{W}}_{41}-{\rho }_3{\mathrm{D}}_1({\mathbb{W}}_3,{\mathbb{W}}_4)-{\omega }_3{\mathbb{W}}_{41}\right),\hfill \\ & & ⋮=⋮,\hfill \\ & & \mathcal{L}{\mathbb{W}}_{1(\upsilon +1)}\left(t\right)=\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left(\gamma {\mathbb{W}}_{1\left(\upsilon \right)}-\frac{\gamma {\left({\mathrm{A}}_{\left(\upsilon \right)}\left({\mathbb{W}}_1\right)\right)}^2}{\mathcal{L}}-{\rho }_1{\mathrm{B}}_{\left(\upsilon \right)}({\mathbb{W}}_1,{\mathbb{W}}_2)\right),\hfill \\ & & \mathcal{L}{\mathbb{W}}_{2(\upsilon +1)}\left(t\right)=\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left({\rho }_1{\mathrm{𝔟}}_1{\mathrm{B}}_{\left(\upsilon \right)}({\mathbb{W}}_1,{\mathbb{W}}_2)-{\rho }_2{\mathrm{C}}_{\left(\upsilon \right)}({\mathbb{W}}_2,{\mathbb{W}}_3)-{\omega }_1{\mathbb{W}}_{2\left(\upsilon \right)}\right),\hfill \\ & & \mathcal{L}{\mathbb{W}}_{3(\upsilon +1)}\left(t\right)=\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left({\rho }_2{\mathrm{𝔟}}_2{\mathrm{C}}_{\left(\upsilon \right)}({\mathbb{W}}_2,{\mathbb{W}}_3)+{\rho }_3{\mathrm{𝔟}}_3{\mathrm{D}}_{\left(\upsilon \right)}({\mathbb{W}}_3,{\mathbb{W}}_4)+\wp {\mathbb{W}}_{4\left(\upsilon \right)}-{\omega }_2{\mathbb{W}}_{3\left(\upsilon \right)}\right),\hfill \\ & & \mathcal{L}{\mathbb{W}}_{4(\upsilon +1)}\left(t\right)=\frac{2(s+\mathbf{𝔡}(1-s\left)\right)}{sa(2-\mathbf{𝔡})}\mathcal{L}\left(\alpha {\mathbb{W}}_{3\left(\upsilon \right)}-\wp {\mathbb{W}}_{4\left(\upsilon \right)}-{\rho }_3{\mathrm{D}}_{\left(\upsilon \right)}({\mathbb{W}}_3,{\mathbb{W}}_4)-{\omega }_3{\mathbb{W}}_{4\left(\upsilon \right)}\right),\upsilon \ge 0.\hfill \end{array}$$

(22)

Now, applying ${\mathcal{L}}^{-1}$, on both sides of (22) and using the following notations for the sake of simplicity,

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathsf{H}}_1=\gamma {\mathbb{W}}_{10\upsilon }-\frac{\gamma {\left({\mathrm{A}}_0\left({\mathbb{W}}_1\right)\right)}^2}{\mathcal{L}}-{\rho }_1{\mathrm{B}}_0({\mathbb{W}}_1,{\mathbb{W}}_2),\hfill \\ & {\mathsf{H}}_2={\rho }_1{\mathrm{𝔟}}_1{\mathrm{B}}_0({\mathbb{W}}_1,{\mathbb{W}}_2)-{\rho }_2{\mathrm{C}}_0({\mathbb{W}}_2,{\mathbb{W}}_3)-{\omega }_1{\mathbb{W}}_{20},\hfill \\ & {\mathsf{H}}_3={\rho }_2{\mathrm{𝔟}}_2{\mathrm{C}}_0({\mathbb{W}}_2,{\mathbb{W}}_3)+{\rho }_3{\mathrm{𝔟}}_3{\mathrm{D}}_0({\mathbb{W}}_3,{\mathbb{W}}_4)+\wp {\mathbb{W}}_{40}-{\omega }_2{\mathbb{W}}_{30},\hfill \\ & {\mathsf{H}}_4=\alpha {\mathbb{W}}_{30}-\wp {\mathbb{W}}_{40}-{\rho }_3{\mathrm{D}}_0({\mathbb{W}}_3,{\mathbb{W}}_4)-{\omega }_3{\mathbb{W}}_{40}.\hfill \end{array}\end{array}$$

(23)

Thus, after using the above values, one can obtain

$$\begin{array}{ccc}& & {\mathbb{W}}_{10}\left(t\right)={\mathbb{W}}_1\left(0\right),{\mathbb{W}}_{20}\left(t\right)={\mathbb{W}}_2\left(0\right),{\mathbb{W}}_{30}\left(t\right)={\mathbb{W}}_3\left(0\right),{\mathbb{W}}_{40}\left(t\right)={\mathbb{W}}_4\left(0\right),\hfill \\ & & {\mathbb{W}}_{11}\left(t\right)=\frac{2(1-\mathbf{𝔡}+\mathbf{𝔡}t)}{2-\mathbf{𝔡}}{\mathfrak{H}}_1,{\mathbb{W}}_{21}\left(t\right)=\frac{2(1-\mathbf{𝔡}+\mathbf{𝔡}t)}{2-\mathbf{𝔡}}{\mathfrak{H}}_2,{\mathbb{W}}_{31}\left(t\right)=\frac{2(1-\mathbf{𝔡}+\mathbf{𝔡}t)}{2-\mathbf{𝔡}}{\mathfrak{H}}_3,\hfill \\ & & {\mathbb{W}}_{41}\left(t\right)=\frac{2(1-\mathbf{𝔡}+\mathbf{𝔡}t)}{2-\mathbf{𝔡}}{\mathfrak{H}}_4;\hfill \end{array}$$

(24)

similarly, the other terms can be calculated, and the final series solution can be written as

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathbb{W}}_{1\upsilon }\left(t\right)={\mathbb{W}}_{10}\left(t\right)+{\mathbb{W}}_{11}\left(t\right)+{\mathbb{W}}_{12}\left(t\right)+\dots ,\hfill \\ & {\mathbb{W}}_{2\upsilon }\left(t\right)={\mathbb{W}}_{20}\left(t\right)+{\mathbb{W}}_{21}\left(t\right)+{\mathbb{W}}_{22}\left(t\right)+\dots ,\hfill \\ & {\mathbb{W}}_{3\upsilon }\left(t\right)={\mathbb{W}}_{30}\left(t\right)+{\mathbb{W}}_{31}\left(t\right)+{\mathbb{W}}_{32}\left(t\right)+\dots ,\hfill \\ & {\mathbb{W}}_{4\upsilon }\left(t\right)={\mathbb{W}}_{40}\left(t\right)+{\mathbb{W}}_{41}\left(t\right)+{\mathbb{W}}_{42}\left(t\right)+\dots .\hfill \end{array}\end{array}$$

(25)

7. Results and Discussion

This section aims to present the numerical simulations of the approximate solution of the system (2). The phase projections are present to analyse the behavior and effects of important parameters as well as fractional order on the dynamics of model (2). The technique used for the analytical approximation is highly converging and can easily be implemented. For the simulation purpose, the initial values are considered as $[{\mathbb{W}}_1,{\mathbb{W}}_2{\mathbb{W}}_3,{\mathbb{W}}_4]=[1,1,1,1]$. The parameters’ values are considered as presented in

Table 1. The parameters and their description of the model (2).

Parameter	Value	Parameter	Value
$\gamma$	1	${\rho }_1$	1
${\rho }_2$	0.25	${\rho }_3$	0.1
${\mathrm{𝔟}}_1$	0.65	${\mathrm{𝔟}}_2$	0.5
${\omega }_1$	0.01	${\omega }_2$	0.2
${\omega }_3$	0.01	℘	0.15
$\mathcal{L}$	100	$\alpha$	0.15
${\mathrm{𝔟}}_3$	0.5

, which are taken from [32]. In Figure 2, Figure 3 and Figure 4, the fractional orders are considered as $(\mathrm{𝔟}lue,1.00),(green,0.99),(red,0.98),(purple,0.97)$. The MATLAB version R2018a has been used for the simulations. The time $t=100$ is considered for the numerical simulations with $dt=0.001$.

In Figure 2a–d, the dynamics of the system state variables ${\mathbb{W}}_1,{\mathbb{W}}_2,{\mathbb{W}}_3\mathrm{and}{\mathbb{W}}_4$ are presented, respectively, vs time t. It can be observed that, at an integer higher fractional order, there are a number of oscillations in the system, and the system behaves chaotically, while, at lower fractional orders, the oscillation amplitudes decreases, resulting in the faster converging towards the equilibrium point. The system state variables becomes stable at $t=38$ with $d=0.97$. Furthermore, in Figure 3, the 2D behaviors of different state variables are demonstrated. In Figure 3a, the dynamics of prey and intermediate populations are presented; similarly, in Figure 3b, the behavior of prey and immature populations is projected while Figure 3c shows the population dynamics of intermediate and immature predators. Moreover, in Figure 4, the 3D dynamics of different state variables are presented in which Figure 4a shows the population behavior of intermediate, mature, and immature predators. Figure 4b,c demonstrates the dynamics of prey, intermediate predator, mature predators and prey, intermediate predator, and immature predators, respectively. From these, the complex nature of the system can be observed.

In Figure 5, the dynamical evolution of the different populations with different initial conditions is presented. The blue color shows the behavior of different populations with initial conditions $[{\mathbb{W}}_1,{\mathbb{W}}_2,{\mathbb{W}}_3,{\mathbb{W}}_4]=10,8,6,4$, green color shows the dynamics with initial conditions $8,6,4,2$, and red curves show the evolution with initial conditions $6,4,6,8$. In Figure 5a, the population dynamics of intermediate, mature and prey, intermediate predators are depicted, respectively, with $\mathrm{𝔡}=1.00$. Similarly, Figure 5b,c shows the population behavior of intermediate, mature and prey, intermediate predators, respectively, with $\mathrm{𝔡}=0.97$. It can be observed that, at lower fractional orders, the complexity in the dynamics is reduced, and the system moves towards stability quickly as compared to higher fractional orders. From the simulations, it is observed that, at higher fractional orders, there are a number of oscillations in the system, and it behaves chaotically, while, at lower fractional orders, the amplitudes of oscillations decrease, resulting in faster converging towards the equilibrium point.

In Figure 6, the dynamical behavior of prey population, mature, and immature predators, respectively. In Figure 6a,b, the blue, green, and red colors show the behavior of different populations with varying the parameter $\mathcal{L}$. From Figure 6c,d, the blue, green, and red colors show the behavior of different populations with varying the parameter $\gamma$. According to the simulations, the system exhibits a variety of oscillations and exhibits chaotic behavior at higher fractional orders, whereas, at lower fractional orders, the oscillations’ amplitudes diminish, and the system rapidly approaches equilibrium.

8. Conclusions

This article investigated the dynamics of the food-web model with three population species in the sense of modified CF derivatives. The operator used for the analysis of the model is the modified version of the CF operator. With the help of fixed point theory, existence and uniqueness solution of the fractional food web model are presented. Analytically, the suggested model is approximated using the Laplace–Adomian technique. The technique used for the analytical approximation is highly converging and can easily be implemented. The solution of the system is presented graphically with different fractional orders as well as different initial conditions. With the help of a graphical representation, it was demonstrated how the system’s parameters and derivative order both had a significant impact. The new modified fractional derivative will be applied to other disease models in the future with memory or hereditary properties.