On Numerical Analysis of Bio-Ethanol Production Model with the Effect of Recycling and Death Rates under Fractal Fractional Operators with Three Different Kernels

Abstract

The main metabolism of yeasts produces bioethanol. Bioethanol, which is produced from biomass and bioenergy crops, has been promoted as one of the most viable alternatives to fossil fuels. The following reaction represents all of the knowledge we have regarding intracellular reactions and their regulatory mechanisms: $biomass+substrates\to ethanol+biomass\left(morecells\right).$ Atangana has suggested new operators based on a combination of fractional and fractal calculus. Fractal-fractional operators (FFOs) have frequently been utilized to investigate the dynamics of a physical problem. In this paper, FFOs are used to investigate a nonlinear mathematical model for ethanol production with three different kernels. Famous fixed point results are employed to show the existence and uniqueness of the solution of the FFO ethanol model under the Mittag–Leffler kernel. The concept of non-linear analysis is utilized to demonstrate the model’s Ulam–Hyres stability. The Adams—Bashforth numerical technique, which is based on the Lagrangian interpolation method, is utilized to find the solution of the model under fractal-fractional operators with three different kernels. The numerical results are simulated with MATLAB-17 for several sets of fractional orders and fractal dimensions to show the relationship between components of ethanol production under new operators in various senses.

1. Introduction

In recent years, scientists and engineers have focused on the development of green, sustainable bioenergies. The creation of a large spectrum of compounds from sustainable derivations, such as biomass or energy, is implied by this idea [1,2]. Processes that convert biomass into usable biomaterial have the potential to enhance the economic value of currently wasted raw materials, while simultaneously lowering the amount of wastewater produced by various businesses. Bioethanol is an essential renewable fuel that can aid in mitigating the negative environmental effects of worldwide fossil fuel consumption. To reduce the amount of greenhouse gases discharged into the environment, ethanol can be added to gasoline as a transportation fuel [3]. Bioethanol is a renewable liquid biofuel that may be used to replace oil-based fossil fuels. Agricultural crops are currently the most popular bioethanol synthesis substrates (sorghum, corn, sugarcane, wheat). Recognizing how a bioreactor’s processes produce bioethanol requires determining the rate at which biomass grows. There are a variety of models that look at the evolution of cell mass development in a bioreactor, and the bulk of them provide formulas for the particular growth rate of cell mass [4]. Many researchers have been curious about the long-term behavior of enlarged models with a continuous flow reactor [5,6]. Ajbar examined the creation of complicated dynamic behaviour of competition in microbes [7]. Cornelli et al. looked at the capacity of yeasts to generate sugars in a range of products manufactured by famous firms in 2013 [8]. Bhowmik et al. [9] developed a model for bioethanol production in 2018. The model of Cornelli et al. was expanded by the author to include the recycling ratio and decay rate.

Here, we consider the bioethanol production model given in [9] by using novel operators. The model consists of three contents: The symbol $\mathbf{S}$ represents the evolution of the substrate, $\mathbf{B}$ represents the biomass, and $\mathbf{E}$ denotes the ethanol. The dimensional bioethanol model is as follows:

$$\left\{\begin{array}{c}\mathcal{V}\frac{d}{dt}\mathbf{S}=\mathcal{F}({\mathbf{S}}_0-\mathbf{S})-\frac{{ϵ}_{max}}{{\mathrm{Y}}_{\frac{\mathbf{X}}{\mathbf{S}}}}{\mathcal{M}}_2(\mathbf{S}, \mathbf{E})\mathbf{X}\mathcal{V},\hfill \\ \mathcal{V}\frac{d}{dt}\mathbf{B}=-\mathcal{F}\mathbf{B}+{ϵ}_{max}{\mathcal{M}}_2(\mathbf{S}, \mathbf{E})\mathbf{B}\mathcal{V}+\mathrm{R}\mathcal{F}(\mathrm{C}-1)\mathbf{B}-{\mathrm{b}}_{\mathrm{H}}\mathbf{B}\mathcal{V},\hfill \\ \mathcal{V}\frac{d}{dt}\mathbf{E}=-\mathcal{F}\mathbf{E}+{\mathrm{Y}}_{\frac{\mathbf{E}}{\mathbf{B}}}{ϵ}_{max}{\mathcal{M}}_2(\mathbf{S}, \mathbf{E})\mathbf{B}\mathcal{V}+\mathrm{R}\mathcal{F}(\mathrm{C}-1)\mathbf{E}+\mathcal{F}\gamma \mathbf{B},\hfill \end{array}\right.$$

(1)

where ${\mathcal{M}}_2$ represents the modified Andrew expression with inhibition by ethanol. The following equations denote the ${\mathcal{M}}_2$ and residence time $\tau$:

$$\begin{array}{cc}\hfill {\mathcal{M}}_2(\mathbf{S}, \mathbf{E})& =\frac{\mathbf{S}}{{\mathcal{K}}_{\mathbf{S}}+\mathbf{S}+{\mathcal{K}}_{\mathbf{E}}{\mathbf{E}}^2},\hfill \\ \hfill \tau & =\frac{\mathcal{V}}{\mathcal{F}}.\hfill \end{array}$$

Variables $\mathbf{S}, \mathbf{B}$, and $\mathbf{E}$ denote the substrate, biomass, and ethanol, respectively. The three components of the investigated model all have non-negative initial values, that is, $\mathbf{S}\left(0\right)={\mathbf{S}}_0\ge 0, \mathbf{B}\left(0\right)={\mathbf{B}}_0\ge 0,$ and $\mathbf{E}\left(0\right)={\mathbf{E}}_0\ge 0.$ The dimensionless version of the model (1) is established by expressing the dimensionless variables as:

$${\mathbf{S}}^{*}=\frac{\mathbf{S}}{{\mathcal{K}}_{\mathbf{S}}}, {\mathbf{B}}^{*}=\frac{\mathbf{B}}{{\mathrm{Y}}_{\frac{\mathbf{E}}{\mathbf{B}}}{\mathcal{K}}_{\mathbf{S}}}, {\mathbf{E}}^{*}=\frac{\mathbf{E}}{{\mathcal{K}}_{\mathbf{S}}}, \mathit{t}={ϵ}_{max}t.$$

Now model (1) in dimensionless form is given by:

$$\left\{\begin{array}{c}\frac{d}{dt}{\mathbf{S}}^{*}\left(\mathit{t}\right)=\frac{{\mathbf{S}}_0^{*}-{\mathbf{S}}^{*}}{{\tau }^{*}}-\frac{{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}},\hfill \\ \frac{d}{dt}{\mathbf{B}}^{*}\left(\mathit{t}\right)=\frac{-{\mathbf{B}}^{*}}{{\tau }^{*}}+\frac{{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}}-{\mathrm{b}}_{\mathrm{H}}^{*}{\mathbf{B}}^{*}+\frac{{\mathrm{R}}^{*}{\mathbf{B}}^{*}}{{\tau }^{*}},\hfill \\ \frac{d}{dt}{\mathbf{E}}^{*}\left(\mathit{t}\right)=\frac{-{\mathbf{E}}^{*}}{{\tau }^{*}}+\frac{{\mathrm{Y}}_2{\mathbf{B}}^{*}}{{\tau }^{*}}+\frac{{\mathrm{Y}}_3{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}}+\frac{{\mathrm{R}}^{*}{\mathbf{E}}^{*}}{{\tau }^{*}}.\hfill \end{array}\right.$$

(2)

All parameters of (2) are non-negative, and ${\mathbf{S}}^{*}\left(0\right)={\mathbf{S}}_0^{*}\ge 0, {\mathbf{B}}^{*}\left(0\right)={\mathbf{B}}_0^{*}\ge 0,$ and ${\mathbf{E}}^{*}\left(0\right)={\mathbf{E}}_0^{*}\ge 0.$ The following is a list of the parameters utilized in the aforementioned systems:

• The symbol $\mathcal{F}$ represents the flow rate by means of reactor.
• ${\mathrm{b}}_{\mathrm{H}}$ represents the coefficient of death.
• ${\mathrm{b}}_{\mathrm{H}}^{*}$ is the dimensionless death rate.
• ${\mathcal{K}}_{\mathbf{S}}$ and ${\mathcal{K}}_{\mathbf{E}}$ are the saturation constant and the inhibition constant by ethanol, respectively.
• ${\mathbf{B}}^{*}$, ${\mathbf{S}}^{*}$, and ${\mathbf{E}}^{*}$ represent the dimensionless biomass, substrate, and ethanol concentrations, respectively.
• The volume of the reactor is denoted by $\mathcal{V}$.
• Time and dimensionless time are denoted by t and $t^{*}$, respectively.
• ${\mathrm{Y}}_{\frac{\mathbf{B}}{\mathbf{S}}}$ and ${\mathrm{Y}}_{\frac{\mathbf{E}}{\mathbf{B}}}$ represent the biomass and ethanol/biomass yield coefficients, respectively.
• $\mathrm{Y}$ denotes the ethanol production’s kinetic constant.
• The rate of specific growth rate is denoted by ${\mathcal{M}}_2(\mathbf{S}, \mathbf{E})$.
• ${\mu }_{max}$ is the maximum specific growth rate.
• The residence time is denoted by $\tau$, while the dimensionless residence time is represented by ${\tau }^{*}$.
• $\mathrm{R}$ is the recycling ratio, which is calculated using volume flow rates.
• Effective recycling is represented by ${\mathrm{R}}^{*}$.

Fractional calculus (FC) has been implemented in different applied disciplines with great success. There are many operators in the literature of FC, but three main operators have been widely used for modeling by researchers. The first is the Caputo operator, the second is the Caputo–Fabrizio (CF) differential operator, and the third is the Atangana–Baleanu (AB) differential operator. These operators are based on power law, exponential decay law, and Mittag–Leffler function, respectively. Mathematical modeling of various phenomena utilizing these operators is highly influenced by the study of fractional differential equations (FDEs) [10,11,12,13,14,15,16,17,18]. Models incorporating FO integral and DEs are shown to be superior to traditional models. Similarly, if the order of the fractal approaches 1, ordinary order differentiation is extended to the idea of fractal differentiation equivalent to the classical derivative. We provide applications of a fractal operator in [19,20]. Recently, Atangana [21] has defined fractal-fractional operators (FFOs) of the function convolution based on three known kernels: the power law, exponential decay law, and the extended Mittag–Leffler function. They take into consideration fractal impact, memory, and non-locality, while combining the concepts of fractal and fractional derivative. The FFOs allow us to understand crossover behavior issues, forgotten memory, and self-similar power-law. FFOs anticipate the complicated behavior of physical issues that are not described by classical operators or FOs [22,23]. Ahmad et al. investigated the complex dynamics of a multi-scroll chaotic system with a Caputo FFO [24]. The authors analyzed a tumor immune model under AB FFO [25]. Owolabi et al. studied a reaction diffusion model under FFO [26]. Akgül et al. investigated magnetohydrodynamics couple stress fluid under FFOs [27]. Some other applications are included in [28,29,30].

Motivated by the previous studies, we examine the ethanol model (2), utilizing the newly constructed fractal-fractional derivatives with three distinct kernels. Using the famous fixed point results of Leray Schauder’s alternative and Banach, we study the existence and uniqueness of the model’s solutions. Additionally, Ulam’s type stability is investigated using nonlinear analysis. Lastly, we implement the Adams–Basforth method to obtain the approximate solutions to the FFO model of ethanol production using three alternative kernels. A graphical representation of the obtained results are depicted to show the dynamics of the components of ethanol production. The following is a summary of the paper’s structure: We start with some fundamental results and definitions in Section 2. Section 3 describes the essential features of the proposed model, such as theoretical results relating to the existence of a solution, and UH stability of the solution. Section 4 is devoted to the Adams–Basforth approach for solving the proposed model. The graphics are explained in detail in Section 5. The conclusion of the paper is provided in Section 6.

2. Preliminaries

Let $\delta$ and $\vartheta$ denote the fractional order and fractal dimension, respectively, and $0\le \delta , \vartheta \le 1$. Let $\mathcal{G}\left(\mathit{t}\right)$ be a continuous and fractal differentiable function on $\left(\mathrm{m}, \mathrm{n}\right).$ We give definitions of FFOs below.

Definition 1

([20]).The Caputo fractal fractional derivative of $\mathcal{G}\left(\mathit{t}\right)$ is defined as:

$${}^{\mathcal{C}}{\mathcal{D}}_{0,t}^{\delta ,\vartheta }\left(\mathcal{G}\left(\mathit{t}\right)\right)=\frac{1}{\Gamma (m-\delta )}{\int }_0^{\mathit{t}}{(\mathit{t}-\lambda )}^{m-\delta -1}\frac{d}{d{\lambda }^{\vartheta }}\mathcal{G}\left(\lambda \right)d\lambda ,$$

where $\frac{d}{d{\lambda }^{\vartheta }}\mathcal{G}\left(\lambda \right)={lim}_{t\to \lambda }\frac{\mathcal{G}\left(\mathit{t}\right)-\mathcal{G}\left(\lambda \right)}{t^{\vartheta }-{\lambda }^{\vartheta }}$.

Definition 2

([20]).The Caputo–Fabrizio fractal fractional derivative of $\mathcal{G}\left(\mathit{t}\right)$ is defined as:

$${}^{\mathcal{C}F}{\mathcal{D}}_{0,t}^{\delta ,\vartheta }\left(\mathcal{G}\left(\mathit{t}\right)\right)=\frac{\mathcal{M}\left(\delta \right)}{1-\delta }{\int }_0^{\mathit{t}}exp\left[-\frac{\delta }{1-\delta }{(\mathit{t}-\lambda )}^{w-\delta -1}\right]\frac{d}{d{\lambda }^{\vartheta }}\mathcal{G}\left(\lambda \right)d\lambda ,$$

where $0<\delta , \vartheta \le w\in \mathbb{N}$ and $\mathcal{M}\left(0\right)=\mathcal{M}\left(1\right)=1.$

Definition 3

([20]).The Atangana–Baleanu fractal fractional derivative of $\mathcal{G}\left(\mathit{t}\right)$ is defined as:

$${}^{\mathcal{ABC}}{\mathcal{D}}_{0,t}^{\delta ,\vartheta }\left(\mathcal{G}\left(\mathit{t}\right)\right)=\frac{\mathcal{AB}\left(\delta \right)}{1-\delta }{\int }_0^{\mathit{t}}{E}_{\delta }\left[-\frac{\delta }{1-\delta }{(\mathit{t}-\lambda )}^{\delta }\right]\frac{d}{d{\lambda }^{\vartheta }}\mathcal{G}\left(\lambda \right)d\lambda ,$$

where $0<\delta , \vartheta \le 1$, and $\mathcal{AB}\left(\delta \right)=1-\delta +\frac{\delta }{\Gamma \left(\delta \right)}.$

Definition 4

([20]).The fractal fractional integral of $\mathcal{G}\left(\mathit{t}\right)$ with power law kernel is given as:

$${}^{\mathcal{C}}{\mathcal{I}}_{0,t}^{\delta ,\vartheta }\mathcal{G}\left(\mathit{t}\right)=\frac{1}{\Gamma \delta }{\int }_0^{\mathit{t}}{(\mathit{t}-\lambda )}^{\delta -1}{\lambda }^{1-\vartheta }\mathcal{G}\left(\lambda \right)d\lambda .$$

Definition 5

([20]).The fractal fractional integral of $\mathcal{G}\left(\mathit{t}\right)$ with exponential decay kernel is defined as:

$${}^{\mathcal{CF}}{\mathcal{I}}_{0,t}^{\delta ,\vartheta }\mathcal{G}\left(\mathit{t}\right)=\frac{\vartheta (1-\delta )t^{\vartheta -1}\mathcal{G}\left(\mathit{t}\right)}{\mathcal{M}\left(\delta \right)}+\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}{\int }_0^{\mathit{t}}{\lambda }^{\delta -1}\mathcal{G}\left(\lambda \right)d\lambda .$$

Definition 6

([20]).The fractal fractional integral of $\mathcal{G}\left(\mathit{t}\right)$ with Mittag–Leffler kernel is defined as:

$${}^{\mathcal{ABC}}{\mathcal{I}}_{0,t}^{\delta ,\vartheta }\mathcal{G}\left(\mathit{t}\right)=\frac{\vartheta (1-\delta )t^{\vartheta -1}\mathcal{G}\left(\mathit{t}\right)}{\mathcal{AB}\left(\delta \right)}+\frac{\delta \vartheta }{\mathcal{AB}\left(\delta \right)}{\int }_0^{\mathit{t}}{\lambda }^{\delta -1}{(\mathit{t}-\lambda )}^{\delta -1}\mathcal{G}\left(\lambda \right)d\lambda .$$

3. Theoretical Results

In this part, we deduce both existence and uniqueness results, as well as the HU-stability results, of the system (2) under FFO in the ABC sense. One may deduce the same results for the Caputo and Caputo–Fabrizio fractal-fractional model of (2). Consider the following:

$$\left\{\begin{array}{c}{}^{\mathcal{ABC}}{\mathcal{D}}_{0,t}^{\delta ,\vartheta }{\mathbf{S}}^{*}\left(\mathit{t}\right)=\frac{{\mathbf{S}}_0^{*}-{\mathbf{S}}^{*}}{{\tau }^{*}}-\frac{{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}},\hfill \\ {}^{\mathcal{ABC}}{\mathcal{D}}_{0,t}^{\delta ,\vartheta }{\mathbf{B}}^{*}\left(\mathit{t}\right)=\frac{-{\mathbf{B}}^{*}}{{\tau }^{*}}+\frac{{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}}-{\mathrm{b}}_{\mathrm{H}}^{*}{\mathbf{B}}^{*}+\frac{{\mathrm{R}}^{*}{\mathbf{B}}^{*}}{{\tau }^{*}},\hfill \\ {}^{\mathcal{ABC}}{\mathcal{D}}_{0,t}^{\delta ,\vartheta }{\mathbf{E}}^{*}\left(\mathit{t}\right)=\frac{-{\mathbf{E}}^{*}}{{\tau }^{*}}+\frac{{\mathrm{Y}}_2{\mathbf{B}}^{*}}{{\tau }^{*}}+\frac{{\mathrm{Y}}_3{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}}+\frac{{\mathrm{R}}^{*}{\mathbf{E}}^{*}}{{\tau }^{*}}.\hfill \end{array}\right.$$

3.1. Existence and Uniqueness Results

With the aid of fixed point theorems, we show that the model under examination has at least one and only one solution. We may formulate the proposed model as follows because the integral is differentiable:

$$\left\{\begin{array}{c}{}_0^{ABR}{D}_{\mathit{t}}^{\delta }{\mathbf{S}}^{*}\left(\mathit{t}\right)=\vartheta {\mathit{t}}^{\vartheta -1}\mathcal{L}(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*}),\hfill \\ {}_0^{ABR}{D}_{\mathit{t}}^{\delta }{\mathbf{B}}^{*}\left(\mathit{t}\right)=\vartheta {\mathit{t}}^{\vartheta -1}\mathcal{W}(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*}),\hfill \\ {}_0^{ABR}{D}_{\mathit{t}}^{\delta }{\mathbf{E}}^{*}\left(\mathit{t}\right)=\vartheta {\mathit{t}}^{\vartheta -1}\mathcal{Z}(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*}),\hfill \end{array}\right.$$

(3)

where

$$\left\{\begin{array}{c}\mathcal{L}(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})=\frac{{\mathbf{S}}_0^{*}-{\mathbf{S}}^{*}}{{\tau }^{*}}-\frac{{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}},\hfill \\ \mathcal{W}(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})=\frac{-{\mathbf{B}}^{*}}{{\tau }^{*}}+\frac{{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}}-{\mathrm{b}}_{\mathrm{H}}^{*}{\mathbf{B}}^{*}+\frac{{\mathrm{R}}^{*}{\mathbf{B}}^{*}}{{\tau }^{*}},\hfill \\ \mathcal{Z}(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})=\frac{-{\mathbf{E}}^{*}}{{\tau }^{*}}+\frac{{\mathrm{Y}}_2{\mathbf{B}}^{*}}{{\tau }^{*}}+\frac{{\mathrm{Y}}_3{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}}+\frac{{\mathrm{R}}^{*}{\mathbf{E}}^{*}}{{\tau }^{*}}.\hfill \end{array}\right.$$

The considered system (3) can be written as

$$\left\{\begin{array}{c}{}_0^{ABR}{D}_{\mathit{t}}^{\delta }\mathcal{O}\left(\mathit{t}\right)=\vartheta {\mathit{t}}^{\vartheta -1}\Lambda \left(t, \mathcal{O}\left(\mathit{t}\right)\right),\hfill \\ \mathcal{O}\left(0\right)={\mathcal{O}}_0.\hfill \end{array}\right.$$

(4)

By replacing ${}_0^{ABR}{D}_{\mathit{t}}^{\delta ,\vartheta }$ with ${}_0^{ABC}{D}_{\mathit{t}}^{\delta ,\vartheta }$ and applying the fractional integral, we get the following:

$$\mathcal{O}\left(\mathit{t}\right)=\mathcal{O}\left(0\right)+\frac{\vartheta t^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}\Lambda \left(t,\mathcal{O}\left(\mathit{t}\right)\right)+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\underset{0}{\overset{\mathit{t}}{\int }}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\vartheta -1}\Lambda \left(\lambda ,\mathcal{O}\left(\lambda \right)\right)d\lambda ,$$

(5)

where

$$\mathcal{O}\left(\mathit{t}\right)=\left\{\begin{array}{c}{\mathbf{S}}^{*}\left(\mathit{t}\right)\hfill \\ {\mathbf{B}}^{*}\left(\mathit{t}\right)\hfill \\ {\mathbf{E}}^{*}\left(\mathit{t}\right)\hfill \end{array}\right.,\mathcal{O}\left(0\right)=\left\{\begin{array}{c}{\mathbf{S}}^{*}\left(0\right)\hfill \\ {\mathbf{B}}^{*}\left(0\right)\hfill \\ {\mathbf{E}}^{*}\left(0\right)\hfill \end{array}\right.,\Lambda \left(t,\mathcal{O}\left(\mathit{t}\right)\right)=\left\{\begin{array}{c}\mathcal{L}(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})\hfill \\ \mathcal{W}(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})\hfill \\ \mathcal{Z}(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})\hfill \end{array}\right..$$

To proceed further, a Banach space is defined as: $\mathbb{B}=\mathsf{C}\times \mathsf{C}\times \mathsf{C},$ where $\mathsf{C}=\mathbf{Q}\left[0,{\mathbb{T}}^{*}\right]$. The norm on the Banach space is defined as

$$∥\mathcal{O}∥=\underset{\mathit{t}\in \left[0,{\mathbb{T}}^{*}\right]}{max}\left|{\mathbf{S}}^{*}\left(\mathit{t}\right)+{\mathbf{B}}^{*}\left(\mathit{t}\right)+{\mathbf{E}}^{*}\left(\mathit{t}\right)\right|.$$

Equation (5) must have a solution if it has a fixed point. Define an operator $\aleph :\mathbb{B}\to \mathbb{B}$. Using (5), one can define ℵ as

$$\aleph \left(\mathcal{O}\right)\left(\mathit{t}\right)=\mathcal{O}\left(0\right)+\frac{\vartheta {\mathit{t}}^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}\Lambda \left(\mathit{t},\mathcal{O}\left(\mathit{t}\right)\right)+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\underset{0}{\overset{\mathit{t}}{\int }}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\vartheta -1}\Lambda \left(\lambda ,\mathcal{O}\left(\lambda \right)\right)d\lambda .$$

(6)

We shift the examined model to a fixed point problem (i.e., $\mathcal{O}=\aleph \left(\mathcal{O}\right)).$ This will help to build a fixed point theory. The following results guarantee the existence of a solution of the model (3).

Theorem 1.

If any $\mathcal{O}\in \mathbb{B},$∃ constants ${\mathcal{T}}_{\Lambda }>0$ and $M_{\Lambda }$ such that

$$\left|\Lambda \left(\mathit{t},\mathcal{O}\left(\mathit{t}\right)\right)\right|\le {\mathcal{T}}_{\Lambda }\left|\mathcal{O}\left(\mathit{t}\right)\right|+M_{\Lambda },$$

(7)

and $\Lambda : \left[0,{\mathbb{T}}^{*}\right]\times \mathbb{B}\to \mathbb{R}$ is a continuous function, then the examined model has at least one solution.

Proof. 

First, we want to prove that ℵ, which is given by (6), is completely continuous. Let $\mathbb{L}=\left\{\mathcal{O}\in \mathbb{B}: ∥\mathcal{O}∥\le \mathcal{R},\mathcal{R}>0\right\}$. For every $\mathcal{O}\in \mathbb{B}$, we have

$$\begin{array}{cc}\hfill ∥\aleph \left(\mathcal{O}\right)∥ =& \underset{0\le \mathit{t}\le {\mathbb{T}}^{*}}{max}\left|\mathcal{O}\left(0\right)+\frac{\vartheta t^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}\Lambda \left(\mathit{t},\mathcal{O}\left(\mathit{t}\right)\right)\right.\hfill \\ \hfill & \left.+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\underset{0}{\overset{\mathit{t}}{\int }}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\vartheta -1}\Lambda \left(\lambda ,\mathcal{O}\left(\lambda \right)\right)d\lambda \right|\hfill \\ \hfill \le & \mathcal{O}\left(0\right)+\frac{\vartheta {\mathbb{T}}^{*}{}^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}\left({\mathcal{T}}_{\Lambda }∥\mathcal{O}∥+M_{\Lambda }\right)\hfill \\ \hfill & +\underset{\mathit{t}\in \left[0,{\mathbb{T}}^{*}\right]}{max}\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\underset{0}{\overset{\mathit{t}}{\int }}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\vartheta -1}\left|\Lambda \left(\lambda ,\mathcal{O}\left(\lambda \right)\right)\right|d\lambda \hfill \\ \hfill \le & \mathcal{O}\left(0\right)+\frac{\vartheta {\mathbb{T}}^{*}{}^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}\left({\mathcal{T}}_{\Lambda }∥\mathcal{O}∥+M_{\Lambda }\right)\hfill \\ \hfill & +\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\left({\mathcal{T}}_{\Lambda }∥\mathcal{O}∥+M_{\Lambda }\right){\mathbb{T}}^{*}{}^{\delta +\vartheta -1}\mathcal{U}(\delta ,\vartheta )\hfill \\ \hfill \le & \mathcal{R},\hfill \end{array}$$

where $\mathcal{U}(\delta ,\vartheta )$ represents the beta function. Thus, we proved the uniform boundedness of the operator ℵ. Now, we have to show the equi-continuity of ℵ. For this, let ${\mathit{t}}_1<{\mathit{t}}_2\le {\mathbb{T}}^{*}$, and we have

$$\begin{array}{cc}\hfill \left|\aleph \left(\mathcal{O}\right)\left({\mathit{t}}_2\right)-\aleph \left(\mathcal{O}\right)\left({\mathit{t}}_1\right)\right|=& \left|\frac{\vartheta {\mathit{t}}_2^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}\Lambda \left({\mathit{t}}_2,\mathcal{O}\left({\mathit{t}}_2\right)\right)+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\underset{0}{\overset{{\mathit{t}}_2}{\int }}{\lambda }^{\vartheta -1}{({\mathit{t}}_2-\lambda )}^{\vartheta -1}\Lambda \left(\lambda ,\mathcal{O}\left(\lambda \right)\right)d\lambda \right.\hfill \\ \hfill & \left.-\frac{\vartheta {\mathit{t}}_1^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}\Lambda \left({\mathit{t}}_1,\mathcal{O}\left({\mathit{t}}_1\right)\right)+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\underset{0}{\overset{{\mathit{t}}_1}{\int }}{\lambda }^{\vartheta -1}{({\mathit{t}}_1-\lambda )}^{\vartheta -1}\Lambda \left(\lambda ,\mathcal{O}\left(\lambda \right)\right)d\lambda \right|\hfill \\ \hfill \le & \frac{\vartheta {\mathit{t}}_2(1-\delta )}{AB\left(\delta \right)}\left({\mathcal{T}}_{\Lambda }\left|\mathcal{O}\left(\mathit{t}\right)\right|+M_{\Lambda }\right)+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\left({\mathcal{T}}_{\Lambda }\left|\mathcal{O}\left(\mathit{t}\right)\right|+M_{\Lambda }\right){\mathit{t}}_2\mathcal{U}(\delta ,\vartheta )\hfill \\ \hfill & -\frac{\vartheta {\mathit{t}}_1(1-\delta )}{AB\left(\delta \right)}\left({\mathcal{T}}_{\Lambda }\left|\mathcal{O}\left(\mathit{t}\right)\right|+M_{\Lambda }\right)-\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\left({\mathcal{T}}_{\Lambda }\left|\mathcal{O}\left(\mathit{t}\right)\right|+M_{\Lambda }\right){\mathit{t}}_1\mathcal{U}(\delta ,\vartheta ).\hfill \end{array}$$

Thus, from the last inequality, it follows that when $\left|\aleph \left(\mathcal{O}\right)\left({\mathit{t}}_2\right)-\aleph \left(\mathcal{O}\right)\left({\mathit{t}}_1\right)\right|\to 0$ as ${\mathit{t}}_1\to {\mathit{t}}_2$. This implies the following:

$$∥\aleph \left(\mathcal{O}\right)\left({\mathit{t}}_2\right)-\aleph \left(\mathcal{O}\right)\left({\mathit{t}}_1\right)∥\to 0,\mathrm{as}{\mathit{t}}_1\to {\mathit{t}}_2.$$

This shows the equi-continuity of ℵ. Using the Arzela–Ascoli theorem, the complete continuity of ℵ follows. As a result, the proposed ethanol model possesses at least one fixed point, according to Schauder’s fixed point result. It follows that at least one solution of the model exists. □

Theorem 2.

If ∀ $\mathcal{O},$$\overline{\mathcal{O}}\in \mathbb{B},$∃ a constant ${\mathcal{A}}_{\Lambda }>0$ such that

$$\left|\Lambda \left(\mathit{t},\mathcal{O}\left(\mathit{t}\right)\right)-\Lambda \left(\mathit{t},\overline{\mathcal{O}}\left(\mathit{t}\right)\right)\right|\le {\mathcal{A}}_{\Lambda }\left|\mathcal{O}\left(\mathit{t}\right)-\overline{\mathcal{O}}\left(\mathit{t}\right)\right|,$$

(8)

and $\rho <1,$ where

$$\rho =\left(\frac{\vartheta \mathbb{T}{}^{*}{}^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}{\mathbb{T}}^{*}{}^{\delta +\vartheta -1}\mathcal{U}(\delta ,\vartheta )\right){\mathcal{A}}_{\Lambda },$$

then the proposed ethanol model has a unique solution.

Proof. 

∀ $\mathcal{O},\overline{\mathcal{O}}\in \mathbb{B},$ and one gets

$$\begin{array}{cc}\hfill ∥\aleph \left(\mathcal{O}\right)-\aleph \left(\overline{\mathcal{O}}\right)∥ =& \underset{\mathit{t}\in \left[0,{\mathbb{T}}^{*}\right]}{max}\left|\frac{\vartheta t^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}\left(\Lambda \left(t,\mathcal{O}\left(\mathit{t}\right)\right)-\Lambda \left(t,\overline{\mathcal{O}}\left(\mathit{t}\right)\right)\right)\right.\hfill \\ \hfill & \left.+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\underset{0}{\overset{\mathit{t}}{\int }}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\vartheta -1}d\lambda \left[\Lambda \left(\lambda ,\mathcal{O}\left(\lambda \right)\right)-\Lambda \left(\lambda ,\overline{\mathcal{O}}\left(\lambda \right)\right)\right]\right|\hfill \\ \hfill \le & \left[\frac{\vartheta {\mathbb{T}}^{*}{}^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}{\mathbb{T}}^{*}{}^{\delta +\vartheta -1}\mathcal{U}(\delta ,\vartheta )\right]∥\mathcal{O}-\overline{\mathcal{O}}∥\hfill \\ \hfill \le & \rho ∥\mathcal{O}-\overline{\mathcal{O}}∥.\hfill \end{array}$$

Since $\rho <1,$ this implies that ℵ is a contraction. According to the Banach fixed result, ℵ has a unique fixed point. As a result, the ethanol system under consideration has a unique solution. □

3.2. Ulam–Hyres Stability

In this section, we will show that the suggested ethanol model is UH-stable.

Definition 7.

The model (3) is UH stable if for any $ϵ>0$ ∃ ${\aleph }_{\delta ,\vartheta }\ge 0$, and for every solution $\mathcal{O}\in \mathbf{Q}(\left[0,{\mathbb{T}}^{*}\right],\mathbb{R})$ of the proposed ethanol model such that

$$\left|{}_0^{FFM}{D}_t^{\delta ,\vartheta }\mathcal{O}\left(\mathit{t}\right)-\Lambda \left(t,\mathcal{O}\left(\mathit{t}\right)\right)\right|\le ϵ,$$

(9)

and ∃ a unique solution $\Re \in \mathbf{Q}(\left[0,{\mathbb{T}}^{*}\right],\mathbb{R})$ of the ethanol model such that

$$\left|\mathcal{O}\left(\mathit{t}\right)-\Re \left(\mathit{t}\right)\right|\le {\aleph }_{\delta ,\vartheta }ϵ,$$

(10)

where $\mathit{t}\in \left[0,{\mathbb{T}}^{*}\right].$

Remark 1.

Let $\mathcal{O}\in \mathsf{C}$ be a solution of (4) if and only if $\exists \chi \in \mathsf{C}$ such that $\chi \left(0\right)=0,$ with the following properties:

• $\left|\chi \left(\mathit{t}\right)\right|\le ϵ,$ for $ϵ>0.$
• ${}_0^{FFM}{D}_t^{\delta ,\vartheta }\mathcal{O}\left(\mathit{t}\right)=\Lambda \left(t,\mathcal{O}\left(\mathit{t}\right)\right)+\chi \left(\mathit{t}\right).$

Lemma 1.

Let $\mathcal{O}\in \mathsf{C}$ be a solution of the following perturbed problem:

$$\begin{array}{cc}\hfill {}_0^{FFM}{D}_t^{\delta ,\vartheta }\mathcal{O}\left(\mathit{t}\right)& =\Lambda \left(t,\mathcal{O}\left(\mathit{t}\right)\right)+\chi \left(\mathit{t}\right),\hfill \\ \hfill \mathcal{O}\left(0\right)& ={\mathcal{O}}_0,\hfill \end{array}$$

then

$$\begin{array}{c}\hfill \left|\mathcal{O}\left(\mathit{t}\right)-\left(\mathcal{O}\left(0\right)+\frac{\vartheta t^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}\Lambda \left(t,\mathcal{O}\left(\mathit{t}\right)\right)+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\underset{0}{\overset{\mathit{t}}{\int }}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\vartheta -1}\Lambda \left(\lambda ,\mathcal{O}\left(\lambda \right)\right)d\lambda \right)\right|\le {\delta }_{\delta ,\vartheta }^{*}ϵ,\end{array}$$

where ${\delta }_{\delta ,\vartheta }^{*}=\frac{\vartheta {\mathbb{T}}^{*}{}^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}{\mathbb{T}}^{*}{}^{\delta +\vartheta -1}\mathcal{U}(\delta ,\vartheta )$.

Proof. 

We can easily achieve the desired result by using Equation (5). □

Lemma 2.

If (8) and Lemma (1) holds, and $\rho <1,$ then the proposed ethanol system is UH stable.

Proof. 

Let $\Re \in \mathbb{B}$ be a unique solution and $\mathcal{O}\in \mathbb{B}$ be any solution of the considered ethanol model. Then,

$$\begin{array}{c}\hfill \left|\mathcal{O}\left(\mathit{t}\right)-\Re \left(\mathit{t}\right)\right|\\ \hfill =& \left|\mathcal{O}\left(\mathit{t}\right)-\left[\Re \left(0\right)+\frac{\vartheta t^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}\Lambda \left(t,\Re \left(\mathit{t}\right)\right)\right.\right.\hfill \\ \hfill & \left.\left.+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\underset{0}{\overset{\mathit{t}}{\int }}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\vartheta -1}\Lambda \left(\lambda ,\Re \left(\lambda \right)\right)d\lambda \right]\right|\hfill \\ \hfill \le & \left|\mathcal{O}\left(\mathit{t}\right)-\left(\mathcal{O}\left(0\right)+\frac{\vartheta t^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}\Lambda \left(t,\mathcal{O}\left(\mathit{t}\right)\right)\right.\right.\hfill \\ \hfill & \left.\left.+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\underset{0}{\overset{\mathit{t}}{\int }}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\vartheta -1}\Lambda \left(\lambda ,\mathcal{O}\left(\lambda \right)\right)d\lambda \right)\right|\hfill \\ \hfill & +\left|\mathcal{O}\left(0\right)+\frac{\vartheta t^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}\Lambda \left(t,\mathcal{O}\left(\mathit{t}\right)\right)+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\underset{0}{\overset{\mathit{t}}{\int }}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\vartheta -1}\Lambda \left(\lambda ,\mathcal{O}\left(\lambda \right)\right)d\lambda \right|\hfill \\ \hfill & -\left|\Re \left(0\right)+\frac{\vartheta t^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}\Lambda \left(t,\Re \left(\mathit{t}\right)\right)+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}\underset{0}{\overset{\mathit{t}}{\int }}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\vartheta -1}\Lambda \left(\lambda ,\Re \left(\lambda \right)\right)d\lambda \right|\hfill \\ \hfill \le & {\delta }_{\delta ,\vartheta }^{*}ϵ+\left(\frac{\vartheta {\mathbb{T}}^{*}{}^{*}{}^{\vartheta -1}(1-\delta )}{AB\left(\delta \right)}+\frac{\delta \vartheta }{AB\left(\delta \right)\mathrm{Y}\left(\delta \right)}{\mathbb{T}}^{*}{}^{*}{}^{\delta +\vartheta -1}\mathcal{U}(\delta ,\vartheta )\right){\mathcal{A}}_{\Lambda }\left|\mathcal{O}\left(\mathit{t}\right)-\Re \left(\mathit{t}\right)\right|\hfill \\ \hfill \le & {\delta }_{\delta ,\vartheta }^{*}ϵ+\rho \left|\mathcal{O}\left(\mathit{t}\right)-\Re \left(\mathit{t}\right)\right|.\hfill \end{array}$$

As a result, one may state the following:

$$∥\mathcal{O}-\Re ∥ \le {\delta }_{\delta ,\vartheta }^{*}ϵ+\rho ∥\mathcal{O}-\Re ∥$$

The above equation may be written as follows:

$$∥\mathcal{O}-\Re ∥ \le {\aleph }_{\delta ,\vartheta }ϵ$$

where ${\aleph }_{\delta ,\vartheta }=\frac{{\delta }_{\delta ,\vartheta }^{*}}{1-\rho }$. This proves UH stability of the ethanol model. □

4. Numerical Schemes

We use the Adams–Bashforth technique to develop numerical techniques for the considered ethanol model using the Caputo, CF, and AB FFOs.

4.1. Adams–Bashforth Method for Power Law Kernal

Consider the model (2) as:

$$\left\{\begin{array}{c}{}^{\mathcal{C}}{\mathcal{D}}_{0,t}^{\delta ,\vartheta }{\mathbf{S}}^{*}\left(\mathit{t}\right)=\frac{{\mathbf{S}}_0^{*}-{\mathbf{S}}^{*}}{{\tau }^{*}}-\frac{{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}},\hfill \\ {}^{\mathcal{C}}{\mathcal{D}}_{0,t}^{\delta ,\vartheta }{\mathbf{B}}^{*}\left(\mathit{t}\right)=\frac{-{\mathbf{B}}^{*}}{{\tau }^{*}}+\frac{{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}}-{\mathrm{b}}_{\mathrm{H}}^{*}{\mathbf{B}}^{*}+\frac{{\mathrm{R}}^{*}{\mathbf{B}}^{*}}{{\tau }^{*}},\hfill \\ {}^{\mathcal{C}}{\mathcal{D}}_{0,t}^{\delta ,\vartheta }{\mathbf{E}}^{*}\left(\mathit{t}\right)=\frac{-{\mathbf{E}}^{*}}{{\tau }^{*}}+\frac{{\mathrm{Y}}_2{\mathbf{B}}^{*}}{{\tau }^{*}}+\frac{{\mathrm{Y}}_3{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}}+\frac{{\mathrm{R}}^{*}{\mathbf{E}}^{*}}{{\tau }^{*}}.\hfill \end{array}\right.$$

(11)

As the fractional integral is differentiable, we may write the given system as:

$$\left\{\begin{array}{c}{}_0^{RL}{D}_{\mathit{t}}^{\delta }{\mathbf{S}}^{*}\left(\mathit{t}\right)=\vartheta {\mathit{t}}^{\vartheta -1}\mathcal{L}(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*}),\hfill \\ {}_0^{RL}{D}_{\mathit{t}}^{\delta }{\mathbf{B}}^{*}\left(\mathit{t}\right)=\vartheta {\mathit{t}}^{\vartheta -1}\mathcal{W}(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*}),\hfill \\ {}_0^{RL}{D}_{\mathit{t}}^{\delta }{\mathbf{E}}^{*}\left(\mathit{t}\right)=\vartheta {\mathit{t}}^{\vartheta -1}\mathcal{Z}(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*}).\hfill \end{array}\right.$$

(12)

To utilize the integer-order initial conditions, we substitute the RL derivative with the Caputo derivative. The fractional integral is then applied on both sides, giving us

$$\left\{\begin{array}{c}{\mathbf{S}}^{*}\left(\mathit{t}\right)={\mathbf{S}}^{*}\left(0\right)+\frac{\vartheta }{\Gamma \left(\delta \right)}{\int }_0^{\mathit{t}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{B}}^{*}\left(\mathit{t}\right)={\mathbf{B}}^{*}\left(0\right)+\frac{\vartheta }{\Gamma \left(\delta \right)}{\int }_0^{\mathit{t}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_2(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{E}}^{*}\left(\mathit{t}\right)={\mathbf{E}}^{*}\left(0\right)+\frac{\vartheta }{\Gamma \left(\delta \right)}{\int }_0^{\mathit{t}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_3(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \end{array}\right.$$

(13)

where

$$\left\{\begin{array}{c}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})=\frac{{\mathbf{S}}_0^{*}-{\mathbf{S}}^{*}}{{\tau }^{*}}-\frac{{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}}\hfill \\ {\mathcal{S}}_2(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})=\frac{-{\mathbf{B}}^{*}}{{\tau }^{*}}+\frac{{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}}-{\mathrm{b}}_{\mathrm{H}}^{*}{\mathbf{B}}^{*}+\frac{{\mathrm{R}}^{*}{\mathbf{B}}^{*}}{{\tau }^{*}},\hfill \\ {\mathcal{S}}_3(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})=\frac{-{\mathbf{E}}^{*}}{{\tau }^{*}}+\frac{{\mathrm{Y}}_2{\mathbf{B}}^{*}}{{\tau }^{*}}+\frac{{\mathrm{Y}}_3{\mathbf{S}}^{*}{\mathbf{B}}^{*}}{1+{\mathbf{S}}^{*}+{\mathrm{Y}}_1{\mathbf{E}}^{{*}^2}}+\frac{{\mathrm{R}}^{*}{\mathbf{E}}^{*}}{{\tau }^{*}},\hfill \end{array}\right.$$

(14)

We now derive the numerical results of the given model at $\mathit{t}={\mathit{t}}_{\mathrm{c}+1},$ so (13) becomes

$$\left\{\begin{array}{cc}{\mathbf{S}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{S}}^{*}{}^0+\frac{\vartheta }{\Gamma \left(\delta \right)}{\int }_0^{t_{\mathrm{c}+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{B}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{B}}^{*}{}^0+\frac{\vartheta }{\Gamma \left(\delta \right)}{\int }_0^{t_{\mathrm{c}+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_2(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{E}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{E}}^{*}{}^0+\frac{\vartheta }{\Gamma \left(\delta \right)}{\int }_0^{t_{\mathrm{c}+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_3(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda .\hfill \end{array}\right.$$

(15)

Then, the approximation of the above system is given below:

$$\left\{\begin{array}{cc}{\mathbf{S}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{S}}^{*}{}^0+\frac{\vartheta }{\Gamma \left(\delta \right)}\sum _{s=0}^{\mathrm{c}}{\int }_{t_s}^{t_{s+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{B}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{B}}^{*}{}^0+\frac{\vartheta }{\Gamma \left(\delta \right)}\sum _{s=0}^{\mathrm{c}}{\int }_{t_s}^{t_{s+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_2(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{E}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{E}}^{*}{}^0+\frac{\vartheta }{\Gamma \left(\delta \right)}\sum _{s=0}^{\mathrm{c}}{\int }_{t_s}^{t_{s+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_3(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda .\hfill \end{array}\right.$$

(16)

We use Lagrangian piece-wise interpolation in the interval $\left[t_s,t_{s+1}\right]$ to approximate the kernel inside the integrals as follows:

$$\left\{\begin{array}{cc}{\mathcal{P}}_s\left(\lambda \right)\hfill & =\frac{\lambda -t_{s-1}}{t_s-t_{s-1}}{t}_s^{\vartheta -1}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^s, {\mathbf{B}}^{*}{}^s, {\mathbf{E}}^{*}{}^s)\hfill \\ \hfill & -\frac{\lambda -t_s}{t_s-t_{s-1}}{t}_{s-1}^{\vartheta -1}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^{s-1}, {\mathbf{B}}^{*}{}^{s-1}, {\mathbf{E}}^{*}{}^{s-1}),\hfill \\ {\mathcal{Q}}_s\left(\lambda \right)\hfill & =\frac{\lambda -t_{s-1}}{t_s-t_{s-1}}{t}_s^{\vartheta -1}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^s, {\mathbf{B}}^{*}{}^s, {\mathbf{E}}^{*}{}^s)\hfill \\ \hfill & -\frac{\lambda -t_s}{t_s-t_{s-1}}{t}_{j-1}^{\vartheta -1}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^{s-1}, {\mathbf{B}}^{*}{}^{s-1}, {\mathbf{E}}^{*}{}^{s-1}),\hfill \\ {\mathcal{R}}_s\left(\lambda \right)\hfill & =\frac{\lambda -t_{s-1}}{t_s-t_{s-1}}{t}_j^{\vartheta -1}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^s, {\mathbf{B}}^{*}{}^s, {\mathbf{E}}^{*}{}^s)\hfill \\ \hfill & -\frac{\lambda -t_s}{t_s-t_{s-1}}{t}_{j-1}^{\vartheta -1}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^{s-1}, {\mathbf{B}}^{*}{}^{s-1}, {\mathbf{E}}^{*}{}^{s-1}).\hfill \end{array}\right.$$

(17)

Thus, the system (16) becomes

$$\left\{\begin{array}{cc}{\mathbf{S}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{S}}^{*}{}^0+\frac{\vartheta }{\Gamma \left(\delta \right)}\sum _{s=0}^{\mathrm{c}}{\int }_{t_s}^{t_{s+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{P}}_s\left(\lambda \right)d\lambda ,\hfill \\ {\mathbf{B}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{B}}^{*}{}^0+\frac{\vartheta }{\Gamma \left(\delta \right)}\sum _{s=0}^{\mathrm{c}}{\int }_{t_s}^{t_{s+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{Q}}_s\left(\lambda \right)d\lambda ,\hfill \\ {\mathbf{E}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{E}}^{*}{}^0+\frac{\vartheta }{\Gamma \left(\delta \right)}\sum _{s=0}^{\mathrm{c}}{\int }_{t_s}^{t_{s+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{R}}_s\left(\lambda \right)d\lambda .\hfill \end{array}\right.$$

(18)

We get the following numerical scheme by solving the integrals of the right-hand sides.

$$\left\{\begin{array}{cc}{\mathbf{S}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{S}}^{*}{}^0+\frac{\vartheta {\left(\mathsf{\Delta }t\right)}^{\delta }}{\Gamma (\delta +2)}\sum _{s=0}^{\mathrm{c}}\left[t_{s-1}^{\vartheta }{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^s, {\mathbf{B}}^{*}{}^s, {\mathbf{E}}^{*}{}^s)\right.\hfill \\ \hfill & \times \left({(1+\mathrm{c}-s)}^{\delta }(2+\mathrm{c}+\delta -s)-{(\mathrm{c}-s)}^{\delta }(\mathrm{c}+2\delta +2-s)\right)\hfill \\ \hfill & -t_{s-1}^{\vartheta -1}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^{s-1}, {\mathbf{B}}^{*}{}^{s-1}, {\mathbf{E}}^{*}{}^{s-1})\hfill \\ \hfill & \left.\times \left({(1+\mathrm{c}-s)}^{\delta +1}-{(\mathrm{c}-s)}^{\delta }(\mathrm{c}-s+1+\delta )\right)\right],\hfill \\ {\mathbf{B}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{B}}^{*}{}^0+\frac{\vartheta {\left(\mathsf{\Delta }t\right)}^{\delta }}{\Gamma (\delta +2)}\sum _{s=0}^{\mathrm{c}}\left[t_{s-1}^{\vartheta }{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^s, {\mathbf{B}}^{*}{}^s, {\mathbf{E}}^{*}{}^s)\right.\hfill \\ \hfill & \times \left({(1+\mathrm{c}-s)}^{\delta }(2+\mathrm{c}+\delta -s)-{(\mathrm{c}-s)}^{\delta }(\mathrm{c}+2\delta +2-s)\right)\hfill \\ \hfill & -t_{s-1}^{\vartheta -1}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^{s-1}, {\mathbf{B}}^{*}{}^{s-1}, {\mathbf{E}}^{*}{}^{s-1})\hfill \\ \hfill & \left.\times \left({(1+\mathrm{c}-s)}^{\delta +1}-{(\mathrm{c}-s)}^{\delta }(\mathrm{c}-s+1+\delta )\right)\right],\hfill \\ {\mathbf{E}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{E}}^{*}{}^0+\frac{\vartheta {\left(\mathsf{\Delta }t\right)}^{\delta }}{\Gamma (\delta +2)}\sum _{s=0}^{\mathrm{c}}\left[t_{s-1}^{\vartheta }{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^s, {\mathbf{B}}^{*}{}^s, {\mathbf{E}}^{*}{}^s)\right.\hfill \\ \hfill & \times \left({(1+\mathrm{c}-s)}^{\delta }(2+\mathrm{c}+\delta -s)-{(\mathrm{c}-s)}^{\delta }(\mathrm{c}+2\delta +2-s)\right)\hfill \\ \hfill & -t_{s-1}^{\vartheta -1}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^{s-1}, {\mathbf{B}}^{*}{}^{s-1}, {\mathbf{E}}^{*}{}^{s-1})\hfill \\ \hfill & \left.\times \left({(1+\mathrm{c}-s)}^{\delta +1}-{(\mathrm{c}-s)}^{\delta }(\mathrm{c}-s+1+\delta )\right)\right].\hfill \end{array}\right.$$

(19)

4.2. Adams–Bashforth Method for Exponential-Decay Kernel

Consider the model (2) under CF operator as:

$$\left\{\begin{array}{c}{}^{\mathcal{CF}}{\mathcal{D}}_{0,t}^{\delta }{\mathbf{S}}^{*}\left(\mathit{t}\right)=\vartheta t^{\vartheta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*}),\hfill \\ {}^{\mathcal{CF}}{\mathcal{D}}_{0,t}^{\delta }{\mathbf{B}}^{*}\left(\mathit{t}\right)=\vartheta t^{\vartheta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*}),\hfill \\ {}^{\mathcal{CF}}{\mathcal{D}}_{0,t}^{\delta }{\mathbf{E}}^{*}\left(\mathit{t}\right)=\vartheta t^{\vartheta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*}).\hfill \end{array}\right.$$

(20)

Consequently, we get

$$\left\{\begin{array}{cc}{\mathbf{S}}^{*}\left(\mathit{t}\right)\hfill & ={\mathbf{S}}^{*}\left(0\right)+\frac{\vartheta t^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}{\int }_0^{\mathit{t}}{\lambda }^{\vartheta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{B}}^{*}\left(\mathit{t}\right)\hfill & ={\mathbf{B}}^{*}\left(0\right)+\frac{\vartheta t^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}{\int }_0^{\mathit{t}}{\lambda }^{\vartheta -1}{\mathcal{S}}_2(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{E}}^{*}\left(\mathit{t}\right)\hfill & ={\mathbf{E}}^{*}\left(0\right)+\frac{\vartheta t^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})\hfill \\ \hfill & \frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}{\int }_0^{\mathit{t}}{\lambda }^{\vartheta -1}{\mathcal{S}}_3(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda .\hfill \end{array}\right.$$

(21)

Now, we deduce the numerical technique at $t=t_{\mathrm{c}+1}$. Therefore,

$$\left\{\begin{array}{cc}{\mathbf{S}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{S}}^{*}{}^0+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}{\int }_0^{t_{\mathrm{c}+1}}{\lambda }^{\vartheta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{B}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{B}}^{*}{}^0+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}{\int }_0^{t_{\mathrm{c}+1}}{\lambda }^{\vartheta -1}{\mathcal{S}}_2(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{E}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{E}}^{*}{}^0+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & \frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}{\int }_0^{t_{\mathrm{c}+1}}{\lambda }^{\vartheta -1}{\mathcal{S}}_3(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda .\hfill \end{array}\right.$$

(22)

Taking the difference between the consecutive terms, we get

$$\left\{\begin{array}{cc}{\mathbf{S}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{S}}^{*}{}^{\mathrm{c}}+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & -\frac{\vartheta {\mathit{t}}_{\mathrm{c}-1}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}-1}, {\mathbf{B}}^{*}{}^{\mathrm{c}-1}, {\mathbf{E}}^{*}{}^{\mathrm{c}-1})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}{\int }_{t_{\mathrm{c}}}^{t_{\mathrm{c}+1}}{\lambda }^{\vartheta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{B}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{B}}^{*}{}^{\mathrm{c}}+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & -\frac{\vartheta {\mathit{t}}_{\mathrm{c}-1}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}-1}, {\mathbf{B}}^{*}{}^{\mathrm{c}-1}, {\mathbf{E}}^{*}{}^{\mathrm{c}-1})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}{\int }_{t_{\mathrm{c}}}^{t_{\mathrm{c}+1}}{\lambda }^{\vartheta -1}{\mathcal{S}}_2(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{E}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{E}}^{*}{}^{\mathrm{c}}+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & -\frac{\vartheta {\mathit{t}}_{\mathrm{c}-1}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}-1}, {\mathbf{B}}^{*}{}^{\mathrm{c}-1}, {\mathbf{E}}^{*}{}^{\mathrm{c}-1})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}{\int }_{t_{\mathrm{c}}}^{t_{\mathrm{c}+1}}{\lambda }^{\vartheta -1}{\mathcal{S}}_3(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda .\hfill \end{array}\right.$$

(23)

Now, using integration and Lagrange polynomial interpolation, we get

$$\left\{\begin{array}{cc}{\mathbf{S}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{S}}^{*}{}^{\mathrm{c}}+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & -\frac{\vartheta {\mathit{t}}_{\mathrm{c}-1}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}-1}, {\mathbf{B}}^{*}{}^{\mathrm{c}-1}, {\mathbf{E}}^{*}{}^{\mathrm{c}-1})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}\frac{3}{2}\left(\mathsf{\Delta }t\right){\mathit{t}}_{\mathrm{c}}^{\vartheta -1}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & -\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}\frac{\mathsf{\Delta }t}{2}{\mathit{t}}_{\mathrm{c}-1}^{\vartheta -1}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}-1}, {\mathbf{B}}^{*}{}^{\mathrm{c}-1}, {\mathbf{E}}^{*}{}^{\mathrm{c}-1}),\hfill \\ {\mathbf{B}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{B}}^{*}{}^{\mathrm{c}}+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & -\frac{\vartheta {\mathit{t}}_{\mathrm{c}-1}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}-1}, {\mathbf{B}}^{*}{}^{\mathrm{c}-1}, {\mathbf{E}}^{*}{}^{\mathrm{c}-1})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}\frac{3}{2}\left(\mathsf{\Delta }t\right){\mathit{t}}_{\mathrm{c}}^{\vartheta -1}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & -\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}\frac{\mathsf{\Delta }t}{2}{\mathit{t}}_{\mathrm{c}-1}^{\vartheta -1}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}-1}, {\mathbf{B}}^{*}{}^{\mathrm{c}-1}, {\mathbf{E}}^{*}{}^{\mathrm{c}-1}),\hfill \\ {\mathbf{E}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{E}}^{*}{}^{\mathrm{c}}+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & -\frac{\vartheta {\mathit{t}}_{\mathrm{c}-1}^{\vartheta -1}(1-\delta )}{\mathcal{M}\left(\delta \right)}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}-1}, {\mathbf{B}}^{*}{}^{\mathrm{c}-1}, {\mathbf{E}}^{*}{}^{\mathrm{c}-1})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}\frac{3}{2}\left(\mathsf{\Delta }t\right){\mathit{t}}_{\mathrm{c}}^{\vartheta -1}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & -\frac{\delta \vartheta }{\mathcal{M}\left(\delta \right)}\frac{\mathsf{\Delta }t}{2}{\mathit{t}}_{\mathrm{c}-1}^{\vartheta -1}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}-1}, {\mathbf{B}}^{*}{}^{\mathrm{c}-1}, {\mathbf{E}}^{*}{}^{\mathrm{c}-1}).\hfill \end{array}\right.$$

(24)

4.3. Adams–Bashforth Method for Mittag–Leffler Type Kernel

Here, we consider the system (2) in terms of the AB FFO as:

$$\left\{\begin{array}{c}{}^{\mathcal{ABR}}{\mathcal{D}}_{0,t}^{\delta }{\mathbf{S}}^{*}\left(\mathit{t}\right)=\vartheta t^{\vartheta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*}),\hfill \\ {}^{\mathcal{ABR}}{\mathcal{D}}_{0,t}^{\delta }{\mathbf{B}}^{*}\left(\mathit{t}\right)=\vartheta t^{\vartheta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*}),\hfill \\ {}^{\mathcal{ABR}}{\mathcal{D}}_{0,t}^{\delta }{\mathbf{E}}^{*}\left(\mathit{t}\right)=\vartheta t^{\vartheta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*}).\hfill \end{array}\right.$$

(25)

Taking the given system in the ABC sense and applying the AB integral, we have

$$\left\{\begin{array}{cc}{\mathbf{S}}^{*}\left(\mathit{t}\right)\hfill & ={\mathbf{S}}^{*}\left(0\right)+\frac{\vartheta t^{\vartheta -1}(1-\delta )}{\mathcal{AB}\left(\delta \right)}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{AB}\left(\delta \right)\Gamma \left(\delta \right)}{\int }_0^{\mathit{t}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{B}}^{*}\left(\mathit{t}\right)\hfill & ={\mathbf{B}}^{*}\left(0\right)+\frac{\vartheta t^{\vartheta -1}(1-\delta )}{\mathcal{AB}\left(\delta \right)}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{AB}\left(\delta \right)\Gamma \left(\delta \right)}{\int }_0^{\mathit{t}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_2(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{E}}^{*}\left(\mathit{t}\right)\hfill & ={\mathbf{E}}^{*}\left(0\right)+\frac{\vartheta t^{\vartheta -1}(1-\delta )}{\mathcal{AB}\left(\delta \right)}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{AB}\left(\delta \right)\Gamma \left(\delta \right)}{\int }_0^{\mathit{t}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_3(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda .\hfill \end{array}\right.$$

(26)

Now, at $t=t_{\mathrm{c}+1},$ we get

$$\left\{\begin{array}{cc}{\mathbf{S}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{S}}^{*}{}^0+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{AB}\left(\delta \right)}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{AB}\left(\delta \right)\Gamma \left(\delta \right)}{\int }_0^{t_{\mathrm{c}+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*},{\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{B}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{B}}^{*}{}^0+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{AB}\left(\delta \right)}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{AB}\left(\delta \right)\Gamma \left(\delta \right)}{\int }_0^{t_{\mathrm{c}+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_2(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{E}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{E}}^{*}{}^0+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{AB}\left(\delta \right)}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{AB}\left(\delta \right)\Gamma \left(\delta \right)}{\int }_0^{t_{\mathrm{c}+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_3(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda .\hfill \end{array}\right.,$$

(27)

Using the approximation of the integrals, we get

$$\left\{\begin{array}{cc}{\mathbf{S}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{S}}^{*}{}^0+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{AB}\left(\delta \right)}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{AB}\left(\delta \right)\Gamma \left(\delta \right)}\sum _{s=0}^{\mathrm{c}}{\int }_{t_s}^{t_{s+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_1(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{B}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{B}}^{*}{}^0+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{AB}\left(\delta \right)}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{AB}\left(\delta \right)\Gamma \left(\delta \right)}\sum _{s=0}^{\mathrm{c}}{\int }_{t_s}^{t_{s+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_2(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda ,\hfill \\ {\mathbf{E}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{E}}^{*}{}^0+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{AB}\left(\delta \right)}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})\hfill \\ \hfill & +\frac{\delta \vartheta }{\mathcal{AB}\left(\delta \right)\Gamma \left(\delta \right)}\sum _{s=0}^{\mathrm{c}}{\int }_{t_s}^{t_{s+1}}{\lambda }^{\vartheta -1}{(\mathit{t}-\lambda )}^{\delta -1}{\mathcal{S}}_3(\lambda , {\mathbf{S}}^{*}, {\mathbf{B}}^{*}, {\mathbf{E}}^{*})d\lambda .\hfill \end{array}\right.$$

(28)

Now, utilizing Lagrangian polynomial piece-wise interpolation, one can achieve the following

$$\left\{\begin{array}{cc}{\mathbf{S}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{S}}^{*}{}^0+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{AB}\left(\delta \right)}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})+\frac{\vartheta {\left(\mathsf{\Delta }t\right)}^{\delta }}{\mathcal{AB}\left(\delta \right)\Gamma (\delta +2)}\hfill \\ \hfill & \times \sum _{s=0}^{\mathrm{c}}\left[t_s^{\vartheta -1}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^s, {\mathbf{B}}^{*}{}^s, {\mathbf{E}}^{*}{}^s)\right.\hfill \\ \hfill & \times \left({(1+\mathrm{c}-s)}^{\delta }(2+\mathrm{c}+\delta -s)-{(\mathrm{c}-s)}^{\delta }(2\delta +\mathrm{c}-s+2)\right)\hfill \\ \hfill & -t_{s-1}^{\vartheta -1}{\mathcal{S}}_1(\mathit{t}, {\mathbf{S}}^{*}{}^{s-1}, {\mathbf{B}}^{*}{}^{s-1}, {\mathbf{E}}^{*}{}^{s-1})\hfill \\ \hfill & \left.\times \left({(1+\mathrm{c}-s)}^{\delta +1}-{(\mathrm{c}-s)}^{\delta }(1+\delta +\mathrm{c}-s)\right)\right],\hfill \\ {\mathbf{B}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{B}}^{*}{}^0+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{AB}\left(\delta \right)}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})+\frac{\vartheta {\left(\mathsf{\Delta }t\right)}^{\delta }}{\mathcal{AB}\left(\delta \right)\Gamma (\delta +2)}\hfill \\ \hfill & \times \sum _{s=0}^{\mathrm{c}}\left[t_s^{\vartheta -1}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^s, {\mathbf{B}}^{*}{}^s, {\mathbf{E}}^{*}{}^s)\right.\hfill \\ \hfill & \times \left({(1+\mathrm{c}-s)}^{\delta }(2+\mathrm{c}+\delta -s)-{(\mathrm{c}-s)}^{\delta }(2\delta +\mathrm{c}-s+2)\right)\hfill \\ \hfill & -t_{s-1}^{\vartheta -1}{\mathcal{S}}_2(\mathit{t}, {\mathbf{S}}^{*}{}^{s-1}, {\mathbf{B}}^{*}{}^{s-1}, {\mathbf{E}}^{*}{}^{s-1})\hfill \\ \hfill & \left.\times \left({(1+\mathrm{c}-s)}^{\delta +1}-{(\mathrm{c}-s)}^{\delta }(1+\delta +\mathrm{c}-s)\right)\right],\hfill \\ {\mathbf{E}}^{*}{}^{\mathrm{c}+1}\hfill & ={\mathbf{E}}^{*}{}^0+\frac{\vartheta {\mathit{t}}_{\mathrm{c}}^{\vartheta -1}(1-\delta )}{\mathcal{AB}\left(\delta \right)}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^{\mathrm{c}}, {\mathbf{B}}^{*}{}^{\mathrm{c}}, {\mathbf{E}}^{*}{}^{\mathrm{c}})+\frac{\vartheta {\left(\mathsf{\Delta }t\right)}^{\delta }}{\mathcal{AB}\left(\delta \right)\Gamma (\delta +2)}\hfill \\ \hfill & \times \sum _{s=0}^{\mathrm{c}}\left[t_s^{\vartheta -1}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^s, {\mathbf{B}}^{*}{}^s, {\mathbf{E}}^{*}{}^s)\right.\hfill \\ \hfill & \times \left({(1+\mathrm{c}-s)}^{\delta }(2+\mathrm{c}+\delta -s)-{(\mathrm{c}-s)}^{\delta }(2\delta +\mathrm{c}-s+2)\right)\hfill \\ \hfill & -t_{s-1}^{\vartheta -1}{\mathcal{S}}_3(\mathit{t}, {\mathbf{S}}^{*}{}^{s-1}, {\mathbf{B}}^{*}{}^{s-1}, {\mathbf{E}}^{*}{}^{s-1})\hfill \\ \hfill & \left.\times \left({(1+\mathrm{c}-s)}^{\delta +1}-{(\mathrm{c}-s)}^{\delta }(1+\delta +\mathrm{c}-s)\right)\right].\hfill \end{array}\right.$$

(29)

5. Graphical Illustrations

In this section, we use the rate of change of substrate, biomass, and ethanol to simulate the proposed dimensionless nonlinear model. We look at kinetics with sugar concentrations of 100–250 g/L, which are often used in ethanol production. For the graphical representation, we utilise the parameter values and initial conditions as follows: $\mathrm{Y}=0.338$$\frac{g_{ethanol}}{g_{biomass}}$, ${\mu }_{max}=0.333h^{-1}$, ${\mathrm{b}}_{\mathrm{H}}=0.000916h^{-1}$, ${\mathrm{Y}}_{\frac{\mathbf{E}}{\mathbf{B}}}=3.817g_{ethanol}/g_{biomass}$, ${\mathrm{Y}}_{\frac{\mathbf{B}}{\mathbf{S}}}=0.054g_{biomass}/g_{ethanol}$, ${\mathcal{K}}_{\mathbf{E}}=0.048L{g}_{substrate}/g_{ethanol}^2$, ${\mathcal{K}}_{\mathbf{S}}=0.032g_{substrate}/L$, ${\mathbf{S}}_0^{*}=100$, ${\mathbf{B}}_0^{*}=0.5$, and ${\mathbf{S}}_0^{*}=0.3$. Other parameters in system (2) have the following values: ${\mathrm{b}}_{\mathrm{H}}^{*}=\frac{{\mathrm{b}}_{\mathrm{H}}}{{\mu }_{max}}$, ${\mathrm{Y}}_1={\mathcal{K}}_{\mathbf{E}}{\mathcal{K}}_{\mathbf{S}}$, ${\tau }^{*}=\frac{\mathcal{V}{\mu }_{max}}{\mathcal{F}}$, ${\mathrm{Y}}_2=\mathrm{Y}{\mathrm{Y}}_{\frac{\mathbf{B}}{\mathbf{S}}}$, and ${\mathrm{Y}}_3={\mathrm{Y}}_{\frac{\mathbf{E}}{\mathbf{B}}}{\mathrm{Y}}_{\frac{\mathbf{B}}{\mathbf{S}}}$. We show how the substrate, biomass, and ethanol evolve for various $\delta$ and $\vartheta$ values. Figure 1 and Figure 2 represent the graphs of the substrate, biomass, and ethanol concentration for various $\delta$ and $\vartheta =0.8$and$0.9$, respectively, under Caputo FFO. Figure 3 and Figure 4 show the behavior of the considered model for various values of $\delta$ and $\vartheta =0.8$ and $0.9$ in the sense of CF FFO. Figure 5 and Figure 6 are the graphical representation of the proposed model for different values of $\delta$ and $\vartheta =0.8$ and $0.9$ in the sense of AB FFO. Reactors with no or continuous substrate, biomass, or ethanol are considered. The sudden rise in feed concentration caused by the intake of feed may be observed in the substrate concentration. The bacteria then utilize the substrate, which causes it to flatten out to a constant. The concentrations of biomass and ethanol rise rapidly at the initial stage, then over time, they become stable. We have observed the effect of the recycle parameter on the evolution of the proposed model. The statistics in Figure 7 show that increasing recycling increases the production of biomass and methane. As a result, the substrates necessary for a successful biorefinery are depleted. We can observe that the substrate concentration increases as the recycling ratio rises with a fixed time of t = 100. The influence of the death rate on ethanol production dynamics may then be observed. To investigate the influence of the parameter on the dynamical solutions of the proposed model, we begin by changing the death rate using physically and economically significant values. For various selections of death rate, we see a similar rapid rise and decay of substrate in Figure 8. As the death rate increases, the solution of substrate, biomass, and ethanol concentrations move to a steady state over the time. As illustrated in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6, the influence of $\vartheta$ and $\delta$ on each other is also observed. As $\delta$ or $\vartheta$ increases, the slower the process of the substrate, biomass, and ethanol, and vice versa. Furthermore, changing the kernel of the operator produces a small change in the simulations of the proposed model, as seen in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6. The effect of the various kernels on the substrate, biomass, and ethanol concentrations may be shown in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6. From Figure 1 and Figure 2, it can be seen that the rate of increase or decrease in the concentration of the three components of the proposed model is highly sensitive to the fractional order. A small change in the fractional order produces high variation in the dynamics of the substrate, biomass, and ethanol concentration. The peaks or lowest values of the three compartments of the considered model are quickly achieved at a higher fractional order. In addition, stability occurs at a higher fractional order. In Figure 3 and Figure 4, it is observed that the rate of increase or decrease in the concentration of the three components of the proposed model is less sensitive to the fractional order. Varying the fractional order causes little variation in the evolution of substrate, biomass, and ethanol concentrations. The rate of increasing or decreasing in the substrate, biomass, and ethanol concentrations is the same in Figure 5 and Figure 6, up to the peak or lower value. The variation in the dynamics of the three classes occurs after achieving their maximum or minimum value. Overall, we can conclude that the kernel used in FFOs has a significant impact on the ethanol production process. Regarding the fractal dimension, the considered model has different behavior. The evolution of the components of the model changes with changing fractal dimensions.

6. Conclusions

In this paper, we have used the more generalized operators to investigate the ethanol production model with the impacts of recycling and the death rate. We have employed fractal-fractional operators with three different kernels. The existence and uniqueness of the solution of the suggested ethanol in the AB sense have been derived with help of Leray–Schauder and Banach’s fixed point theorem. We have presented the UH stability of the considered model in the sense of the AB operator. We have obtained numerical schemes through the Adams–Bashforth method for three Caputo, CF, and AB fractal-fractional operators. We have used MATLAB-17 (2017b) to simulate the achieved results in order to visualize the evolution of the ethanol model under different FFOs. Graphs have been used to examine the effect of fractal dimension on fractional order. The dynamics of the model’s many compartments are affected by changes in fractal dimension. Further, we have observed the effect of recycling and death rate on the production of ethanol from the substrate. From the numerical simulations, we have noticed that there is variation in the graphical representation of the considered model when we change the kernel of the FFOs. The numerical simulations also revealed that fractal-fractional orders are able to capture more information than usual fractional derivatives due to the fractal dimension. To give a suggestion and direction, we would suggest that the current model can be developed further by introducing a component that captures the random behavior caused by a major source of uncertainty, which usually propagates in time. When such component is added, then the model obtained will be governed by stochastic FDEs. There are many studies and applications of stochastic FDEs in the modelling of a physical process [31,32,33]. Some other future suggestions and extensions of the current model are as follows:

• A rigorous bifurcation investigation of the steady-state solutions with respect to different parameters.
• Investigation of chaotic behaviour.
• Sensitivity and controllability of the considered model.