Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

114 
https://internationalpubls.com 

Revolutionizing Electrochemical kinetics: Homogeneous Redox Catalysis 

unveiled through an innovative mathematical framework using the 

Homotopy Perturbation Method 
 

Saranya.K1, Dr.R. Angel Joy2 

1Research Scholar, Department of Mathematics, Sri GVG Visalakshi College for women, Udumalpet. 

2Associate Professor & HOD of Mathematics, Sri GVG Visalakshi College for women, Udumalpet. 

 

Article History: 

Received: 25-05-2024 

Revised: 11-07-2024 

Accepted: 24-07-2024 

Abstract: 

Introduction: In Redox homogeneous catalysis the catalyst couple merely 

plays the role of an electron in contrast with chemical catalysis. The 

mathematical modelling of homogeneous redox catalysis of electrochemical 

reactions for the EC scheme is described by the system of second order non-

linear differential equations. This model was initially designed using the 

classical differential equations and it is extended to the Caputo fractional 

derivative (FDE’s) of order 𝛼 in particularly sequential fractional case. The 

aim is to obtain approximate analytical solutions of the system of  FDE’s 

using Homotopy Perturbation method (HPM) and analyze the impact of 

various parameters with different order of 𝛼. In addition, the solution of the 

sequential dynamic equations yields the solution of the corresponding 

integer-order differential equations as a special case. 

Keywords: Electrochemical reaction, Fractional Differential Equations, 

Homotopy Perturbation Method. 

 

 

1. Introduction 

Fractional derivatives provide a special instrument for the description of the memory effects and 

hereditary properties of various materials and processes. In the comparison between fractional 

derivatives and integer-order derivatives, these effects are not taken into account. K.B. Oldham et al. 

[2] mentioned in their book that applications in physics, chemistry, and engineering play an 

outstanding role in the development of the subject of applied fractional calculus. In [1], M. Caputo 

published his book in 1969, in which he systematically used the original definition of fractional 

differentiation for formulating and solving various types of problems in various applications like 

viscoelasticity, rheological properties of rocks, and many other fields. See [14, 16, 21]. They 

observed in their research of sequential fractional differential equations with boundary value 

problems and experiments that the use of half-order derivatives and integrals led to a formulation of 

certain electro-chemical problems that was more useful than the classical approach. And also see 

[22], where the author mentions that we can compute the numerical solution of the sequential Caputo 

fractional equations with boundary conditions tend to the corresponding solution of the integer 

boundary value problems. Otherwise, it can be noted that it is non-sequential. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

115 
https://internationalpubls.com 

In the early 1990s, the perturbation method was one of the most effective methods. However, the 

solutions dealt with small parameters. The development of a new method called the Homotopy 

Perturbation Method (HPM) has been done to overcome this problem. Ji-Huan He was the first to 

introduce the homotopy perturbation method in 1998. He introduced this as a means of resolving 

problems in linear, non-linear, and coupled differential equations in ordinary and partial differential 

equations with boundary or initial conditions. HPM is a semi-analytical method that is widely used 

by mathematicians and engineers. This method is a highly effective and convenient way to overcome 

the difficulties of traditional methods. 

Redox reactions generally involve the transfer of electrons between species. Redox reactions are part 

of everyday activities such as photosynthesis, respiration, coal combustion, and fertilizer production 

and use. Homogeneous redox catalysis involves analyzing and predicting the effectiveness of 

catalysis as a function of the potential separation of the catalyst, substrate, and stationary and quasi-

stationary methods. C.P. Andrieux et al. [3-5] discussed the mathematical model of the homogeneous 

redox catalysis of electrochemical reactions. A mathematical model that includes boundary 

conditions features a system of second-order nonlinear differential equations.               

This paper investigates the dynamics of the Electrochemical mechanism model in classical second-

order non-linear differential equations with boundary value conditions [4] initially and then updated 

it to Caputo Sequential fractional equations. The Homotopy Perturbation method (HPM) is used to 

derive the approximate analytical solutions of the FDE system in fractional order. To analyze the 

effect of various values of each parameter with different order 𝛼  on the concentration profile. 

Additionally, the HPM solutions of the EC mechanisms of fractional order tend to the integer-order 

HPM solution of the EC mechanism as a special case 𝛼 → 1. 

Nomenclature 

𝑃, 𝑄 Catalyst Couple. 

𝐴 Substrate. 

𝐵 Product. 

𝑐𝑃 , 𝑐𝑄 , 𝑐𝐴, 𝑐𝐵 Concentration of  𝑃, 𝑄, 𝐴, 𝐵. 

𝐷𝑝, 𝐷𝑄, 𝐷𝐴, 𝐷𝐵 Diffusion coefficients of 𝑃, 𝑄, 𝐴, 𝐵 in 𝑐𝑚2𝑠−1. 

𝐾, 𝐾1, 𝐾2 Rate  constants. 

𝜆, 𝜆1, 𝜆2  Dimensionless rate parameters. 

𝑡 Time. 

𝑥 Distance to the Electrode. 

𝐸 Electrode Potential. 

𝐸𝑃𝑄
𝑜  Standard potential of the Catalyst. 

𝑆     Electrode Surface area. 

𝛿    Polarography of the Diffusion layer thickness. 

𝑣     Kinetic viscosity in 𝑐𝑚 𝑠−1. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

116 
https://internationalpubls.com 

 

 

 

 

 

 

2.Sequential Fractional Derivatives 

This approach is based on the observation [2, 10] that the 𝑚-th order differentiation is simply a series 

of first order differentiations. That is, 

𝑑𝑚𝜑(𝑡)

𝑑𝑡𝑚
= (

𝑑

𝑑𝑡
 .

𝑑

𝑑𝑡
… (𝑚 𝑡𝑖𝑚𝑒𝑠)) 𝜑(𝑡) 

K.S.Miller and B.Ross first called the generalized fractional differentiation defined by 𝐷𝛼 is the 

Riemann-Lioville fractional derivative and Sequential differentiation. That is a suitable method for 

replacing 
𝑑

𝑑𝑡
 with the derivative of non-integer order 𝐷𝛼 where 0 ≤ 𝛼 ≤ 1. Therefore, the general  

fractional differentiation series is written as 

𝐷𝑚𝛼𝜑(𝑡) = (𝐷𝛼 . 𝐷𝛼 … (𝑚 𝑡𝑖𝑚𝑒𝑠))𝜑(𝑡) 

3. Preliminaries 

In this section, some basic definitions that will be used in the remaining sections are given, 

Definition 1. A real valued function 𝜑(𝑡), defined on [𝑎, 𝑏] lies in the function space ℂ𝜇[𝑎, 𝑏], 𝜇 ∈

ℝ, if there exists a real number 𝑘 > 𝜇 such that 𝜑(𝑡) = (𝑡 − 𝑎)𝑘�̅�(𝑡), with �̅� ∈ ℂ[𝑎, 𝑏], and it is 

said to be in the space ℂ𝜇
𝑛 if and only if 𝜑(𝑛) ∈ ℂ𝜇, 𝑛 is positive integer number with zero. It is clear 

that ℂ𝜇1
⊂ ℂ𝜇2

 for 𝜇1 ≤ 𝜇2.   

Definition 2. The Caputo fractional derivative of order  𝛼 ∈ ℝ is defined as  

𝐷𝑡
𝛼

0
𝐶 𝜑(𝑡) =  

1

𝛤(𝑛 − 𝛼)
∫ (𝑡 − 𝜏)𝑛−𝛼−1𝜑(𝑛)(𝜏)𝑑𝜏

𝑡

0

 

where 𝛼 > 0. If 𝛼 = 𝑛 then 𝐷𝑡
𝛼

0
𝐶 𝜑(𝑡) =

𝑑𝑛𝜑(𝑡)

𝑑𝑡𝑛
. 

Definition 3. Let 𝑚𝛼 > 0, and 𝜑(𝑡) ∶ (0, ∞) → ℝ. Then, the Caputo derivative of 𝜑(𝑡) of order 𝑚𝛼 

is given by  

𝐷𝑡
𝑚𝛼

0
𝐶 𝜑(𝑡) =  

1

𝛤(𝑛 − 𝑚𝛼)
∫ (𝑡 − 𝜏)𝑛−𝑚𝛼−1𝜑(𝑛)(𝜏)𝑑𝜏

𝑡

0

 

where 𝑚 ∈ ℕ such that (𝑛 − 1) < 𝑚𝛼 < 𝑛. In particular, if 𝛼 = 1, then 𝑚𝛼 = 𝑚 is an integer and   

𝐷𝑡
𝑚𝛼

0
𝐶 𝜑(𝑡) = 𝐷𝑡

𝑚
0
𝐶 𝜑(𝑡) and 𝐷𝑡

𝛼
0
𝐶 𝜑(𝑡) = 𝜑1(𝑡).  

𝜔 Angular rotation speed in 𝑟𝑎𝑑 𝑠−1. 

𝜉 Dimensionless potential scale. 

𝜓∞ Dimensionless Plateau current. 

𝑖1 Catalytic Plateau current. 

𝑖𝑑𝑃 Diffusion-controlled Plateau current of 𝑃. 

𝑖𝑑𝐴 Diffusion-controlled Plateau current of 𝐴. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

117 
https://internationalpubls.com 

Definition 4. The Caputo fractional derivative of 𝜑(𝑡) of order 𝑚𝛼 for (𝑛 − 1) < 𝑚𝛼 < 𝑛 is said to 

be the sequential Caputo fractional derivative of order 𝛼 if the relation is written as 

𝐷𝑡
𝑚𝛼

0
𝐶 𝜑(𝑡) = 𝐷𝑡

𝛼
0
𝐶 ( 𝐷𝑡

(𝑚−1)𝛼
0
𝐶 )𝜑(𝑡) 

holds for 𝑚 = 2,3, … 

Definition 5. The Riemann-Lioville fractional integral operator (ℐ𝛼) of order 𝛼 > 0, of a function 

𝜑(𝑡) is defined as  

ℐ𝛼𝜑(𝑡) =
1

𝛤(𝛼)
∫ (𝑡 − 𝜏)𝛼−1𝜑(𝜏)𝑑𝜏

𝑡

0

          (𝛼 > 0), 

ℐ0𝜑(𝑡) = 𝜑(𝑡)                                                                     

where  𝛤(𝛼) is the well-known gamma function. 

Some of the properties for ℐ𝛼 and 𝐷𝛼 can be derived as follows 

1) 𝐷𝑡
𝛼

0
𝐶 𝜑𝛿 =

𝛤(𝛿+1)

𝛤(𝛿−𝛼+1)
𝜑𝛿−𝛼,        𝛼 > 0,   𝛿 > −1, 𝜑 > 0,     

2) ℐ𝛼𝜑𝛿 =  
𝛤(𝛿+1)

𝛤(𝛿+𝛼+1)
𝜑𝛿+𝛼,         𝛼 > 0,   𝛿 > −1, 𝜑 > 0,      

3) ℐ𝛼ℐ𝜗𝜑(𝑡) = ℐ𝜗ℐ𝛼𝜑(𝑡) = ℐ𝛼+𝜗𝜑(𝑡),       𝛼, 𝜗 > 0,   

4) 𝐷𝑏
𝛼

𝑎
𝐶  ℐ𝛼𝜑(𝑡) = 𝜑(𝑡),       𝑎 ≤ 𝑡 ≤ 𝑏,   

5)  ℐ𝛼 𝐷𝑡
𝛼

0
𝐶 𝜑(𝑡) = 𝜑(𝑡) − ∑ 𝜑(𝑘)(0+)

𝑡𝑘

𝑘!
 ,

𝜇−1
𝑘          𝜇 − 1 < 𝛼 < 𝜇, 𝑡 > 0. 

Lemma 1. If 𝑚 − 1 < 𝛼 ≤ 𝑚 is the order of Caputo fractional derivative 𝐷𝑏
𝛼

𝑎
𝐶  , then it is consistent 

with the integer-order derivative 
𝑑𝑚

𝑑𝑡𝑚 for 𝑚 ∈ ℕ. 

Proof. For 𝜑 ∈ ℂ𝑚+1([0, ∞]),then 

                  𝐷𝑡
𝛼 𝜑0

𝐶 (𝑡) =
1

𝛤(𝑚 − 𝛼)
∫

  𝜑(𝑚)(𝜏)

(𝑡 − 𝜏)1−(𝑚−𝛼)
𝑑𝜏

𝑡

0

 

                                     =
1

𝛤(𝑚 − 𝛼)
[−

(𝑡 − 𝜏)𝑚−𝛼

𝑚 − 𝛼
 𝜑(𝑚)(𝜏)|𝜏=0

𝜏=𝑡 + ∫
(𝑡 − 𝜏)𝑚−𝛼

𝑚 − 𝛼
 𝜑(𝑚+1)(𝜏)𝑑𝜏

𝑡

0

] 

                                     =
1

𝛤(𝑚 − 𝛼 + 1)
[0 + 𝑡𝑚−𝛼𝜑(𝑚)(0) + ∫ (𝑡 − 𝜏)𝑚−𝛼 𝜑(𝑚+1)(𝜏)𝑑𝜏

𝑡

0

] 

and taking the limit: 𝛼 ∈ ℝ → 𝑚 ∈ ℕ with 𝑚 − 1 ≤ 𝛼 ≤ 𝑚  yields: 

                               =
1

𝛤(1)
[𝜑(𝑚)(0) + ∫ 𝜑(𝑚+1)(𝜏)𝑑𝜏

𝑡

0

] = 𝜑(𝑚)(0) + 𝜑(𝑚)(𝑡) − 𝜑(𝑚)(0) =
𝑑𝑚𝜑

𝑑𝑡𝑚
(𝑡). 

4. Mathematical formulation of the problem 

The simplest of the reaction sequences in which chemical reactions for Electrochemical (EC) 

mechanism can be represented in [4] 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

118 
https://internationalpubls.com 

Electrode process 

 𝐴 + 1𝑒 ⇌ 𝐵                                                                                                                                (1) 

Redox catalyzed process 

 𝑃 + 1𝑒 ⇌ 𝑄(𝐸𝑃𝑄
0 )                                                                                                                    (2) 

            𝑄 + 𝐴  
𝐾1

⇌
𝐾2

  𝑃 + 𝐵                                                                                                                        (3) 

              𝐵
𝐾
→ 𝐶                                                                                                                            (4) 

A reduction process of the oxidized form 𝑃 of the catalyst couple 𝑃 𝑄⁄  is introduced into the 

solution. Here, 𝑃 𝑄⁄  is introduced to fulfill the following conditions: the standard potential 𝐸𝑃𝑄
𝑜   is 

positive for the reduction potential of the substrate 𝐴; electron alteration between the electrode and  

𝑃 or 𝑄 is fast, both 𝑃 and 𝑄 are chemically stable. 

In [3, 4] C.P. Andrieux et.al. treat the system as governed by a set of differential equations with 

boundary conditions and a stationary-state assumption regarding B: 

         𝐾1𝑐𝐴𝑐𝑄 = 𝐾𝑐𝐵 + 𝐾2𝑐𝐵𝑐𝑃                                                                                 (5) 

Introduce the  dimensionless variables as follows:     

𝑦 = 𝑥 𝛿⁄ , 𝑢 =
𝑐𝑃

𝑐𝑃
0  ,        𝑣 =

𝑐𝑄

𝑐𝑃
0  ,         𝑤 =

𝑐𝐴

𝑐𝑃
0  ,        𝑏 =

𝑐𝐵

𝑐𝑃
0  ,       𝛾 =

𝑐𝐴
0

𝑐𝑃
0  , 

𝜉 = −(𝐹 𝑅𝑇⁄ )(𝐸 − 𝐸𝑃𝑄
0 ).                                                                      (6) 

The dimensionless rate parameters[4] 

           𝜆1 =
𝐾1𝑐𝑝

0𝛿2

𝐷
,       𝜆2 =

𝐾2𝑐𝑝
0𝛿2

𝐷
        𝜆 =

𝐾𝛿2

𝐷
                                                                  (7) 

Substituting (6) , (7) in [4] gives 

     
𝑑2𝑢

𝑑𝑦2
= −𝜆1𝑤𝑣 + 𝜆2𝑏𝑢                                                                                                                  (8) 

     
𝑑2𝑣

𝑑𝑦2
= 𝜆1𝑤𝑣 − 𝜆2𝑏𝑢                                                                                                                      (9) 

     
𝑑2𝑤

𝑑𝑦2
= 𝜆1𝑤𝑣 − 𝜆2𝑏𝑢                                                                                                                     (10) 

Using  (6) and (7) in (5)  

     𝜆𝑏 = 𝜆1𝑤𝑣 − 𝜆2𝑏𝑢  or 𝑏 =
𝜆1𝑤𝑣

𝜆+𝜆2𝑢
                                                     (11) 

when  𝑦 = 1:  𝑢 = 1 , 𝑣 = 𝑏 = 0 ,       𝑤 = 𝛾                                                                  (12) 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

119 
https://internationalpubls.com 

when 𝑦 = 0: 

                   (
𝑑𝑢

𝑑𝑦
) + (

𝑑𝑣

𝑑𝑦
) = 0 ,    𝑢 = 𝑣 exp(−𝜉),     (

𝑑𝑤

𝑑𝑦
) = − (

𝑑𝑏

𝑑𝑦
) = 0                             (13) 

From equations (8), (9), (11) and (12), It is clear that  𝑢 + 𝑣 = 1. In particular 𝑢0 + 𝑣0 = 1 and 

therefore,taking (13) into: 

 𝑢0 =
1

1 + exp (𝜉)
, 𝑣0 =

1

1 + exp (−𝜉)
                   

Finally, the system of second order non-linear differential equation and substitute for 𝑏 (11)  in 

(8) − (10) 

   
𝑑2𝑢

𝑑𝑦2
+ 𝜆1𝑤𝑣 − 𝐶𝑢𝑣𝑤 = 0                                                                                                              (14) 

  
𝑑2𝑣

𝑑𝑦2
− 𝜆1𝑤𝑣 + 𝐶𝑢𝑣𝑤 = 0                                                                                                               (15) 

   
𝑑2𝑤

𝑑𝑦2
− 𝜆1𝑤𝑣 + 𝐶𝑢𝑣𝑤 = 0                                                                                                             (16) 

with the boundary conditions, 

                                   𝑦 = 1 ∶       𝑢 = 1 , 𝑣 = 0 , 𝑤 = 𝛾                                                           (17)                                                                                                     

 𝑦 = 0 ∶       𝑢 =
1

[1 + 𝑒𝑥𝑝(𝜉)]
= 𝑠 (say) ,

𝑑𝑢

𝑑𝑦
+

𝑑𝑣

𝑑𝑦
= 0 ,   

𝑑𝑤

𝑑𝑦
= 0                                        (18) 

where   
𝜆1𝜆2

𝜆+𝜆2
= 𝐶(say). 

A fractional model in Caputo Sequential derivative. 

       The fractional framework for the EC mechanism (14) − (16) is given as follows: 

𝐷𝑦
2𝛼

0
𝐶𝑆 𝑢(𝑦) + 𝜆1

𝛼𝑤𝑣 − 𝐶𝛼𝑢𝑣𝑤 = 0                                                                                               (19) 

𝐷𝑦
2𝛼

0
𝐶𝑆 𝑣(𝑦) − 𝜆1

𝛼𝑤𝑣 + 𝐶𝛼𝑢𝑣𝑤 = 0                                                                                               (20) 

𝐷𝑦
2𝛼

0
𝐶𝑆 𝑤(𝑦) − 𝜆1

𝛼𝑤𝑣 + 𝐶𝛼𝑢𝑣𝑤 = 0                                                                                              (21) 

with 𝑏 =
𝜆1

𝛼𝑤𝑣

𝜆𝛼+𝜆2
𝛼𝑢

 .                                                                                                        

                        (22)                                        

where 𝐷𝑦
2𝛼

0
𝐶𝑆 represents Caputo-Sequential derivative with Fractional order 2𝛼. 

5. Analysis of pre-equilibrium and stationary-state assumptions.  

This section discusses the equilibrium points of the EC fractional model and the stationary-state of 

the model (14) − (16) 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

120 
https://internationalpubls.com 

(i) Kinetic control by chemical reaction (4) with electrode process(1) acts as a pre-equilibrium. 

In the context of homogeneous redox catalysis, when the follow-up chemical reaction (4) is slower 

than the backward solution electron transfer (3), the overall reaction kinetics are controlled by the 

follow-up chemical reaction, with the reaction (3) acting as a pre-equilibrium. The proceeding 

system (8) − (10) of three differential equations can be reduced into a single differential equation 

that is 

𝑑2𝑢

𝑑𝑦2
+ (

𝜆𝜆1

𝜆2
) [𝜓(𝑦 − 1) + 1 + 𝛾 − 𝑢]

(1 − 𝑢)

𝑢
= 0                                                                    (23) 

where 𝜓 is the dimensionless plateau current depends upon olnly two parameters 𝛾 and                     
𝜆𝜆1

𝜆2
= (𝑘𝑘1/𝑘2)(𝛿2/𝐷). 

In [4], which corresponds to general cases, some of the observations are represented as follows: 

(i) The system depends on two kinetic parameters and 𝛾, and for high values of 
𝜆𝜆1

𝜆2
, the catalytic 

efficiency reaches unity, and the current-potential curves split into two waves.  

(ii) An increase in δ shifts the catalytic efficiency of the system from low to high.  

(iii) The catalytic efficiency is a function of 𝜆1(7), and an increase in 𝑐0 leads to an increase in 

𝜆1, resulting in the function 
𝜆𝜆1

𝜆2
= (𝑘𝑘1/𝑘2)(𝛿2/𝐷) independent of 𝑐0.   

(ii) The stationary-state assumption regarding B.  

Suppose there is no chemical reaction of the electron-transfer, i.e. when 𝜆 = 0. Then the stationary-

state assumption regarding 𝐵 is no longer valid. 

In [4], (i) The stationary-state assumption is valid when the follow-up chemical reaction is slower 

than the backward solution electron transfer, and the catalytic efficiency is low. (ii) The stationary-

state assumption leads to two limiting situations for extreme values of 𝑘/𝑘2, which involve kinetic 

control of the overall reaction either by the forward reaction (3) or by reaction (4)with (3) as a pre-

equilibrium. (iii) The polarization problem can be numerically solved for these two limiting cases as 

well as for the transition between them in the context of stationary or quasi-stationary 

electrochemical methods. (iv) The stationary-state assumption reduces the number of parameters 

needed to describe the reaction kinetics, making it easier to analyze and understand the kinetics of 

electrochemical reactions under homogeneous redox catalysis. (v) The stationary-state assumption 

has limitations and is only valid for certain ranges of catalytic efficiency and rate constants.  

The previously established results concerning pre-equilibrium and the stationary-state assumption are 

relevant when addressing fractional differential equations with orders ranging from 0 to 1. 

6.Basic idea of Homotopy Perturbation method for system of FDEs 

First to illustrate the basic idea of Homotopy Perturbation method, many applications are modelled 

by systems of Fractional Differential equations which can be written in the form 

𝐷𝛼1𝜑1(𝑡) = 𝑓1(𝑡, 𝜑1, 𝜑2, 𝜑3, … , 𝜑𝑛)                                                                                               (25) 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

121 
https://internationalpubls.com 

𝐷𝛼2𝜑2(𝑡) = 𝑓2(𝑡, 𝜑1, 𝜑2, 𝜑3, … , 𝜑𝑛)                                                                                               (26) 

                

.

.

.
  

𝐷𝛼𝑛𝜑𝑛(𝑡) = 𝑓𝑛(𝑡, 𝜑1, 𝜑2, 𝜑3, … , 𝜑𝑛)                                                                                              (27) 

Subject to the initial condtions: 𝜑𝑛(𝑡) = 𝑐𝑘,   𝑘 = 1,2, … , 𝑛. 

where 𝐷𝛼𝑖 is the fractional derivative of order 𝛼𝑖   (𝑚 − 1 < 𝛼𝑖 ≤ 𝑚) and 𝑓𝑖 are arbitrary linear and 

nonlinar functions.  

construct  the following homotopy: 

 𝐷𝛼𝑖𝜑𝑖 = 𝑝𝑓𝑖(𝑡, 𝜑1, 𝜑2, 𝜑3, … , 𝜑𝑛)                                                                                                    (28) 

where 𝑖 = 1,2, … , 𝑛 and 𝑝 is an embedding parameter which is monotonically increasing from 0 to 1. 

Using the parameter 𝑝, we expand the solution of the system (25) − (27) in the form 

𝜑𝑖(𝑡) = 𝜑𝑖0 + 𝑝𝜑𝑖1 + 𝑝2𝜑𝑖2 + 𝑝3𝜑𝑖3 + ⋯                                                                              (29) 

Substituting (29) into (28) and the resultant equation is in terms with the same powers of 𝑝, 

𝑝𝑜:      𝐷𝛼𝑖𝜑𝑖0 = 0, 

𝑝1:      𝐷𝛼𝑖𝜑𝑖1 = 𝑓𝑖1(𝑡, 𝜑10, 𝜑20, 𝜑30, … , 𝜑𝑛0), 

𝑝2:      𝐷𝛼𝑖𝜑𝑖2 = 𝑓𝑖2(𝑡, 𝜑10, 𝜑20, 𝜑30, … , 𝜑𝑛0,  𝜑11, 𝜑21, 𝜑31, … , 𝜑𝑛1), 

𝑝3:      𝐷𝛼𝑖𝜑𝑖3 = 𝑓𝑖3(𝑡, 𝜑10, 𝜑20, 𝜑30, … , 𝜑𝑛0,  𝜑11, 𝜑21, 𝜑31, … , 𝜑𝑛1, 𝜑11,  𝜑22, 𝜑32, … , 𝜑𝑛2), 

etc., where the functions satisfies the following equation 

𝑓𝑖(𝑡, 𝜑10 + 𝑝𝜑11 + 𝑝2𝜑12 + ⋯ , 𝜑𝑛0 + 𝑝𝜑𝑛1 + 𝑝2𝜑𝑛2 + ⋯ ) 

= 𝑓𝑖1(𝑡, 𝜑10, 𝜑20, 𝜑30, … , 𝜑𝑛0) + 𝑝𝑓𝑖2(𝑡, 𝜑10, 𝜑20, 𝜑30, … , 𝜑𝑛0,  𝜑11, 𝜑21, 𝜑31, … , 𝜑𝑛1)          

+ 𝑝2𝑓𝑖3(𝑡, 𝜑10, 𝜑20, 𝜑30, … , 𝜑𝑛0,  𝜑11, 𝜑21, 𝜑31, … , 𝜑𝑛1, 𝜑11,  𝜑22, 𝜑32, … , 𝜑𝑛2) + ⋯ 

It can be easily solve by applying the ℐ𝛼𝑖  that is the inverse of the operator 𝐷𝛼𝑖. By setting 𝑝 = 1 in 

(28), Hence the components 𝜑𝑖𝑘(𝑘 = 0,1,2, … ) of the HPM solution can be determined. 

7.The HPM of solving system of  FDEs (19)-(21) 

Construct the following homotopy 

𝐷𝑦
2𝛼

0
𝐶𝑆 𝑢(𝑦) = 𝑝(−𝜆1

𝛼𝑤𝑣 + 𝐶𝛼𝑢𝑣𝑤)                                                                                             (30)    

𝐷𝑦
2𝛼

0
𝐶𝑆 𝑣(𝑦) = 𝑝(𝜆1

𝛼𝑤𝑣 − 𝐶𝛼𝑢𝑣𝑤)                                                                                                 (31) 

𝐷𝑦
2𝛼

0
𝐶𝑆 𝑤(𝑦) = 𝑝(𝜆1

𝛼𝑤𝑣 − 𝐶𝛼𝑢𝑣𝑤)                                                                                                (32) 

By perturbation technique, 𝑢, 𝑣 and 𝑤 can be written as a following series in 𝑝 : 

𝑢 = 𝑢0 + 𝑝𝑢1 + 𝑝2𝑢2 + ⋯                                                                                                               (33)    

𝑣 = 𝑣0 + 𝑝𝑣1 + 𝑝2𝑣2 + ⋯                                                                                                               (34)   



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

122 
https://internationalpubls.com 

𝑤 = 𝑤0 + 𝑝𝑤1 + 𝑝2𝑤2 + ⋯                                                                                                            (35)    

Substituting (33) − (35) into (30) − (32) and equating terms of the same powers of  𝑝 yields the 

following three sets of equations 

𝑝0      ∶    𝐷𝑦
2𝛼

0
𝐶𝑆  𝑢0  = 0 𝑝0      ∶      𝐷𝑦

2𝛼
0

𝐶𝑆  𝑣0  = 0 

𝑝1      ∶   𝐷𝑦
2𝛼

0
𝐶𝑆  𝑢1 = −𝜆1

𝛼𝑤0𝑣0 + 𝐶𝛼𝑢0𝑣0𝑤0 𝑝1      ∶      𝐷𝑦
2𝛼

0
𝐶𝑆  𝑣1 = 𝜆1

𝛼𝑤0𝑣0 − 𝐶𝛼𝑢0𝑣0𝑤0 

                      ⋮                        ⋮ 

𝑝0      ∶   𝐷𝑦
2𝛼

0
𝐶𝑆  𝑤0  = 0  

𝑝1      ∶   𝐷𝑦
2𝛼

0
𝐶𝑆  𝑤1 = 𝜆1

𝛼𝑤0𝑣0 − 𝐶𝛼𝑢0𝑣0𝑤0  

                     ⋮  

Applying the operator of  ℐ𝛼 which is inverse operator of 𝐷𝛼 in  to 𝑢0, 𝑣0, 𝑤0, 𝑢1, 𝑣1, 𝑤1, …  

𝑢0(𝑦) = (1 − 𝑠)𝑦 + 𝑠,     𝑣0(𝑦) = (𝑠 − 1)𝑦 + (1 − 𝑠),   𝑤0(𝑦) = 𝛾    

𝑢1(𝑦) = 𝛾(1 − 𝑠)[𝐶𝛼𝑠 − 𝜆1
𝛼]

𝑦2𝛼

𝛤(2𝛼 + 1)
+ 𝛾(𝑠 − 1)[𝐶𝛼(𝑠 − 1) − 𝜆1

𝛼 + 𝐶𝛼𝑠]
𝑦2𝛼+1

𝛤(2𝛼 + 2)
 

                            −𝐶𝛼𝛾(1 − 𝑠)2
𝑦2𝛼+2

𝛤(2𝛼 + 3)
 + {−𝛾(1 − 𝑠)[𝐶𝛼𝑠 − 𝜆1

𝛼]
1

𝛤(2𝛼 + 1)

− 𝛾(𝑠 − 1)[𝐶𝛼(𝑠 − 1) − 𝜆1
𝛼 + 𝐶𝛼𝑠]

1

𝛤(2𝛼 + 2)
+ 𝛾(1 − 𝑠)2

2

𝛤(2𝛼 + 3)
}𝑦𝛼 

 

𝑣1(𝑦) = 𝛾(1 − 𝑠)[𝜆1
𝛼 − 𝐶𝛼𝑠]

𝑦2𝛼

𝛤(2𝛼 + 1)
+ 𝛾(𝑠 − 1)[𝜆1

𝛼 − 𝐶𝛼𝑠 − 𝐶𝛼(𝑠 − 1)]
𝑦2𝛼+1

𝛤(2𝛼 + 2)
 

                           +𝐶𝛼𝛾(1 − 𝑠)2
𝑦2𝛼+2

𝛤(2𝛼 + 3)
 + {−𝛾(1 − 𝑠)[𝐶𝛼𝑠 − 𝜆1

𝛼]
1

𝛤(2𝛼 + 1)

− 𝛾(𝑠 − 1)[𝐶𝛼(𝑠 − 1) − 𝜆1
𝛼 + 𝐶𝛼𝑠]

1

𝛤(2𝛼 + 2)
− 𝛾(1 − 𝑠)2

2

𝛤(2𝛼 + 3)
}𝑦𝛼 

 

𝑤1(𝑦) = 𝛾(1 − 𝑠)[𝜆1
𝛼 − 𝐶𝛼𝑠]

𝑦2𝛼

𝛤(2𝛼 + 1)
+ 𝛾(𝑠 − 1)[𝜆1

𝛼 − 𝐶𝛼𝑠 − 𝐶𝛼(𝑠 − 1)]
𝑦2𝛼+1

𝛤(2𝛼 + 2)
 

                       +𝐶𝛼𝛾(1 − 𝑠)2
𝑦2𝛼+2

𝛤(2𝛼 + 3)
 + {−𝛾(1 − 𝑠)[𝐶𝛼𝑠 − 𝜆1

𝛼]
1

𝛤(2𝛼 + 1)

− 𝛾(𝑠 − 1)[𝐶𝛼(𝑠 − 1) − 𝜆1
𝛼 + 𝐶𝛼𝑠]

1

𝛤(2𝛼 + 2)
− 𝛾(1 − 𝑠)2

2

𝛤(2𝛼 + 3)
}𝑦𝛼 

Letting 𝑝 → 1 gives Therefore, 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

123 
https://internationalpubls.com 

𝑢(𝑦) = (1 − 𝑠)𝑦 + 𝑠 + 𝛾(1 − 𝑠)[𝐶𝛼𝑠 − 𝜆1
𝛼] + 𝛾(𝑠 − 1)[𝐶𝛼(𝑠 − 1) − 𝜆1

𝛼 + 𝐶𝛼𝑠] 

                       
𝑦2𝛼+1

 𝛤(2𝛼 + 2)
− 𝐶𝛼𝛾(1 − 𝑠)2

𝑦2𝛼+2

𝛤(2𝛼 + 3)
 + {−𝛾(1 − 𝑠)[𝐶𝛼𝑠 − 𝜆1

𝛼]
1

𝛤(2𝛼 + 1)

− 𝛾(𝑠 − 1)[𝐶𝛼(𝑠 − 1) − 𝜆1
𝛼 + 𝐶𝛼𝑠]

1

𝛤(2𝛼 + 2)

+ 𝛾(1 − 𝑠)2
2

𝛤(2𝛼 + 3)
}𝑦𝛼                                                                                    (36) 

𝑣(𝑦) = (𝑠 − 1)𝑦 + (1 − 𝑠) + 𝛾(1 − 𝑠)[𝜆1
𝛼 − 𝐶𝛼𝑠]

𝑦2𝛼

𝛤(2𝛼 + 1)

+ 𝛾(𝑠 − 1)[𝜆1
𝛼 − 𝐶𝛼𝑠 − 𝐶𝛼(𝑠 − 1)]

𝑦2𝛼+1

𝛤(2𝛼 + 2)
+ 𝐶𝛼𝛾(1 − 𝑠)2 

𝑦2𝛼+2

𝛤(2𝛼 + 3)
 + {−𝛾(1 − 𝑠)[𝐶𝛼𝑠 − 𝜆1

𝛼]
1

𝛤(2𝛼 + 1)
− 𝛾(𝑠 − 1) 

                      [𝐶𝛼(𝑠 − 1) − 𝜆1
𝛼 + 𝐶𝛼𝑠]

1

𝛤(2𝛼 + 2)
− 𝛾(1 − 𝑠)2

2

𝛤(2𝛼 + 3)
}𝑦𝛼                    (37) 

𝑤(𝑦) = 𝛾 + 𝛾(1 − 𝑠)[𝜆1
𝛼 − 𝐶𝛼𝑠]

𝑦2𝛼

𝛤(2𝛼 + 1)
+ 𝛾(𝑠 − 1)[𝜆1

𝛼 − 𝐶𝛼𝑠 − 𝐶𝛼(𝑠 − 1)]
𝑦2𝛼+1

𝛤(2𝛼 + 2)

+ 𝐶𝛼𝛾(1 − 𝑠)2
𝑦2𝛼+2

𝛤(2𝛼 + 3)
 

                      +{−𝛾(1 − 𝑠)[𝐶𝛼𝑠 − 𝜆1
𝛼]

1

𝛤(2𝛼 + 1)
 − 𝛾(𝑠 − 1)[𝐶𝛼(𝑠 − 1) − 𝜆1

𝛼 + 𝐶𝛼𝑠] 

                       
1

𝛤(2𝛼 + 2)
− 𝛾(1 − 𝑠)2

2

𝛤(2𝛼 + 3)
}                                                                        (38) 

Remark  

Substituting the expressions of 𝑢(𝑦), 𝑣(𝑦) and 𝑤(𝑦) in (22). The concentration of   the product  𝐵 

can be obtained as 𝑏(𝑦).  

8. Analysis and discussion 

    In this section, we established the significance of the parameters 𝜆, 𝐶, 𝛾, 𝜆1, 𝜆2  and  𝜉 on the 

concentration of 𝑃, 𝑄, 𝐴, and 𝐵. 

 

 

 

 

 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

124 
https://internationalpubls.com 

Fig 1 Dimensionless concentration 𝑢 , 𝑣 , 𝑤, 𝑏  versus Dimensionless distance 𝑦  computed using 

(36) to (38), (22) for various values of  𝛼, and for some fixed value of parameters  𝜆 , 𝛾, 𝜆1,  𝜆2 and 

𝜉 . 

   

   

Fig 2 Dimensionless concentration 𝑢 , 𝑣 , 𝑤, 𝑏  versus Dimensionless distance 𝑦  computed using 

(36) to (38), (22) for various values of  𝐶, and for some fixed value of parameters  𝜆 , 𝛾, 𝜆1,  𝜆2 and 

𝜉 . 

   

   



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

125 
https://internationalpubls.com 

Fig 3 Dimensionless concentration 𝑢 , 𝑣 , 𝑤, 𝑏  versus Dimensionless distance 𝑦  computed using 

(36) to (38), (22) for various values of dimensionless potential 𝜉 and for some fixed value of 

parameters  𝜆 , 𝐶, 𝜆1,  𝜆2 and  γ. 

   

   

Fig 4 Dimensionless concentration 𝑢 , 𝑣 , 𝑤, 𝑏  versus Dimensionless distance 𝑦  computed using 

(36) to (38), (22) for various values of 𝜆1, and for some fixed values of parameters 𝜉  , 𝐶,  𝜆2  and  γ. 

   

   



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

126 
https://internationalpubls.com 

Fig 5 The HPM solution of the Caputo sequential system of fractional equations(36) to (38) yields 

the solution of the corresponding integer-order differential equations(𝐴23) , (𝐴24)  and  (𝐴25) as a 

special case 𝛼 → 1. 

   

   

Equations (36) − (38) and the remark represents the analytical expressions for the Dimensionless 

Concentration profile for Catalyst  𝑃 and  𝑄 , and the substrate 𝐴 and the product 𝐵 for the different 

values of the parameters 𝜉, 𝜆1, 𝛾 , 𝐶, 𝜆, 𝜆2 with the different non-integer order. In figures(1), 

demonstrate the Dimensionless concentration of 𝑃, 𝑄, 𝐴, 𝐵 for some fixed values of parameter with 

different values of fractional order 𝛼. This figures gives more accurate results compared to classical 

model. In figures(2), shows that the concentration of the Catalyst  𝑄, the substrate  𝐴 and the product 

𝐵 are directly proportional to 𝐶 but concentration of the Catalyst  𝑃 is inversely proportional to 𝐶. In 

figures(3), Catalyst 𝑃 initially exhibits direct proportionality to potentialξ, changing to inverse 

proportionality at specific points. Similarly, Catalyst 𝑄 initially increases with ξ up to a certain point, 

then decreases. Simultaneously, product 𝐵 undergoes a slight increase followed by a more significant 

decrease. In figures(4), shows that the concentration of the Catalyst  𝑃 is directly proportional to 

𝜆1,but concentration of the Catalyst  𝑄, the substrate 𝐴 and the product 𝐵 are inversely proportional 

to 𝜆1. In figures(5), shows that the fractional HPM yields the solution of the corresponding integer-

order differential equations(𝐴23) , (𝐴24)  and  (𝐴25) as a special case. 

9. Conclusion 

In this paper, a new mathematical model based on the Caputo-Sequential fractional differential 

systems of the EC mechanism was investigated with steady-state assumptions. The fractional model 

of the EC mechanism gives more accurate results as compared to the classical integer-order model. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

127 
https://internationalpubls.com 

The dynamic system of non-linear equations with fractional order has been solved by the semi-

analytical method, that is, HPM, which has very high accuracy. This method is effective and efficient 

not only for classical differential equations but also for fractional derivatives. The influence of the 

various parameters on catalytic efficiency was discussed. Additionally, we demonstrate that the 

sequential system of fractional equations yields the solution of the integer-order differential 

equations as a special case as 𝛼 → 1.    

Appendix A 

The Homotopy Perturbation method of solving system of second order non-linear differential 

equation of EC mechanism 

Here we can solve the system of second order non-linear differential equations (14) − (16) with 

boundary conditions (17) and (18) analytically using Homotopy Perturbation method [4] in the 

following way.  

Using Homotopy for the equations (14) to (16) are 

                        (1 − 𝑝) (
𝑑2𝑢

𝑑𝑦2
) + 𝑝 [

𝑑2𝑢

𝑑𝑦2
+ 𝜆1𝑤𝑣 − 𝐶𝑢𝑣𝑤] = 0                 

                                        

                                                           
𝑑2𝑢

𝑑𝑦2
+ 𝑝𝜆1𝑤𝑣 − 𝑝𝐶𝑢𝑣𝑤 = 0                                                 (𝐴1) 

 

                        (1 − 𝑝) (
𝑑2𝑣

𝑑𝑦2
) + 𝑝 [

𝑑2𝑣

𝑑𝑦2
− 𝜆1𝑤𝑣 + 𝐶𝑢𝑣𝑤] = 0      

                                                           
𝑑2𝑣

𝑑𝑦2
− 𝑝𝜆1𝑤𝑣 + 𝑝𝐶𝑢𝑣𝑤 = 0                                                 (𝐴2) 

                                

                      (1 − 𝑝) (
𝑑2𝑤

𝑑𝑦2
) + 𝑝 [

𝑑2𝑤

𝑑𝑦2
− 𝜆1𝑤𝑣 + 𝐶𝑢𝑣𝑤] = 0   

                                                          
𝑑2𝑤

𝑑𝑦2
− 𝑝𝜆1𝑤𝑣 + 𝑝𝐶𝑢𝑣𝑤 = 0                                                 (𝐴3) 

By Perturbation technique, 𝑢, 𝑣 and 𝑤 can be written as a series in 𝑝 as follows :  

                     𝑢 = 𝑢0 + 𝑝𝑢1 + 𝑝2𝑢2 + ⋯                                     

              (𝐴4)                   

                        𝑣 = 𝑣0 + 𝑝𝑣1 + 𝑝2𝑣2 + ⋯                                 (𝐴5) 

                       𝑤 = 𝑤0 + 𝑝𝑤1 + 𝑝2𝑤2 + ⋯                                             (𝐴6)                         

Equating the identical powers of  𝑝 on both sides of equation  (𝐴1) and (𝐴4) to (𝐴6) gives, 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

128 
https://internationalpubls.com 

𝑝0      ∶          
𝑑2𝑢0

𝑑𝑦2
= 0                                                                                                                         (𝐴7) 

𝑝1      ∶          
𝑑2𝑢1

𝑑𝑦2
+ 𝜆1𝑤0𝑣0 − 𝐶𝑢0𝑤0                                                                                            (𝐴8) 

                      ⋮ 

At 𝑦 = 1, 𝑢0 = 1 ,   𝑢𝑖 = 0   ∀   𝑖 = 1,2,3, …                                                                               (𝐴9) 

At 𝑦 = 0, 𝑢0 = 𝑠 ,   𝑢𝑖 = 0   ∀   𝑖 = 1,2,3, …                                                                             (𝐴10) 

Equating the identical powers of  𝑝 on both sides of equation (𝐴2) and (𝐴4) to (𝐴6) gives, 

𝑝0      ∶          
𝑑2𝑣0

𝑑𝑦2
= 0                                                                                                                      (𝐴11) 

𝑝1      ∶          
𝑑2𝑣1

𝑑𝑦2
− 𝜆1𝑤0𝑣0 + 𝐶𝑢0𝑤0𝑣0 = 0                                                                            (𝐴12) 

                       ⋮ 

At 𝑦 = 1,          𝑣0 = 0    ∀   𝑖 = 0,1,2,3, …                                                                                  (𝐴13) 

At  𝑦 = 0,
𝑑𝑢𝑖

𝑑𝑦
+

𝑑𝑣𝑖

𝑑𝑦
= 0  ∀   𝑖 = 0,1,2,3                                                                            (𝐴14) 

Equating the identical powers of p on both sides of equation (𝐴3) and (𝐴4) to (𝐴6) gives, 

𝑝0      ∶          
𝑑2𝑤0

𝑑𝑦2
= 0                                                                                                                      (𝐴15) 

𝑝1      ∶          
𝑑2𝑤1

𝑑𝑦2
− 𝜆1𝑤0𝑣0 + 𝐶𝑢0𝑤0                                                                                        (𝐴16) 

                       ⋮ 

At 𝑦 = 1,          𝑤0 = 𝛾   𝑤𝑖 = 0  ∀  𝑖 = 1,2,3, …                                                                        (𝐴17) 

At  𝑦 = 0,
𝑑𝑤𝑖

𝑑𝑦
= 0  ∀   𝑖 = 0,1,2,3, …                                                                                 (𝐴18) 

Integrating (𝐴7), (𝐴8), (𝐴11), (𝐴12), (𝐴15), (𝐴16) with respect to y, 

𝑢0(𝑦) = 𝑦(1 − 𝑠) + 𝑠;   𝑣0(𝑦) = (𝑠 − 1)𝑦 + (1 − 𝑠); 𝑤0(𝑦) = 𝛾                       (𝐴19) 

                     

𝑢1(𝑦) =
𝐶𝛾(2𝑠 − 𝑠2 − 1)

12
𝑦4 +

𝛾(2𝐶𝑠2 − 3𝐶𝑠 + 𝐶 − 𝜆1𝑠 + 𝜆1)

6
𝑦3 +

𝛾(𝐶𝑠 − 𝐶𝑠2 − 𝜆1 + 𝜆1𝑠)

2
𝑦2

+
(3𝐶𝑠2𝛾 − 2𝐶𝑠𝛾 − 𝐶𝛾 − 4𝜆1𝑠𝛾 + 4𝜆1𝛾)

12
   

(𝐴20) 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

129 
https://internationalpubls.com 

𝑣1(𝑦) =
𝐶𝛾(𝑠2 − 2𝑠 + 1)

12
𝑦4 +

𝛾(3𝐶𝑠 − 2𝐶𝑠2 − 𝐶 + 𝜆1𝑠 − 𝜆1)

6
𝑦3 +

𝛾(𝐶𝑠2 − 𝐶𝑠 + 𝜆1 − 𝜆1𝑠)

2
𝑦2

+
(2𝐶𝑠𝛾 − 3𝐶𝑠2𝛾 + 𝐶𝛾 − 4𝜆1𝑠𝛾 + 4𝜆1𝛾)

12
 𝑦 

(𝐴21) 

𝑤1(𝑦) =
𝐶𝛾(𝑠2 − 2𝑠 + 1)

12
𝑦4 +

𝛾(3𝐶𝑠 − 2𝐶𝑠2 − 𝐶 + 𝜆1𝑠 − 𝜆1)

6
𝑦3 +

𝛾(𝐶𝑠2 − 𝐶𝑠 + 𝜆1 − 𝜆1𝑠)

2
𝑦2

+
2𝐶𝑠𝛾 − 3𝐶𝑠2𝛾 + 𝐶𝛾 + 4𝜆1𝑠𝛾 − 4𝜆1𝛾

12
 

 (𝐴22) 

Letting 𝑝 → 1 gives, therefore            

𝑢(𝑦) = 𝑠 +
(12 − 12𝑠 + 3𝐶𝑠2𝛾 − 2𝐶𝑠𝛾 − 𝐶𝛾 − 4𝜆1𝑠𝛾 + 4𝜆1𝛾)

12
𝑦 +

𝛾(𝐶𝑠 − 𝐶𝑠2 − 𝜆1 + 𝜆1𝑠)

2
𝑦2

+
𝛾(2𝐶𝑠2 − 3𝐶𝑠 + 𝐶 − 𝜆1𝑠 + 𝜆1)

6
𝑦3 

                              +
𝐶𝛾(2𝑠 − 𝑠2 − 1)

12
 𝑦4                                                                                      (𝐴23) 

𝑣(𝑦) = 1 − 𝑠 −
(12 − 12𝑠 + 3𝐶𝑠2 − 2𝐶𝑠𝛾 − 𝐶𝛾 − 4𝜆1𝑠𝛾 + 4𝜆1𝛾)

12
𝑦 −

𝛾(𝐶𝑠 − 𝐶𝑠2 − 𝜆1 + 𝜆1𝑠)

2
𝑦2

−
𝛾(2𝐶𝑠2 − 3𝐶𝑠 + 𝐶 − 𝜆1𝑠 + 𝜆1)

6
𝑦3 

                                        −
𝐶𝛾(2𝑠 − 𝑠2 − 1)

12
 𝑦4                                                                            (𝐴24) 

𝑤(𝑦) =
(12𝛾 + 2𝐶𝑠𝛾 − 3𝐶𝑠2𝛾 + 𝐶𝛾 + 4𝜆1𝑠𝛾 − 4𝜆1𝛾)

12
−

𝛾(𝐶𝑠 − 𝐶𝑠2 − 𝜆1 + 𝜆1𝑠)

2
𝑦2 

                    −
𝛾(2𝐶𝑠2 − 3𝐶𝑠 + 𝐶 − 𝜆1𝑠 + 𝜆1)

6
𝑦3 −

𝐶𝛾(2𝑠 − 𝑠2 − 1)

12
𝑦4                           (𝐴25) 

The equations (𝐴23) , (𝐴24)  and  (𝐴25)  gives the approximate analytical solution of the system 

(14) − (16). 

Remark  

Substituting the expressions of 𝑢(𝑦), 𝑣(𝑦) and 𝑤(𝑦) in (11). The concentration of  the product  𝐵 

can be obtained as 𝑏(𝑦).   

Reference: 

[1] M.Caputo, ‘Elasticità e Dissipazione’, Zanichelli, Bologna, 1969. 

[2] Oldham.B, Spanier.J, ‘The Fractional Calculus’, Academic Press, New York, NY, USA, London, UK, 1974. 

[3] C.P.Andrieux, J.M.Dumas-Bouchiat, J.M.Saveant, ‘Homogeneous Redox catalysis of electrochemical reaction Part 

I. Introduction’, J. Electroanal. Chem.,87(1978) 39-53. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 6s (2024) 

 

130 
https://internationalpubls.com 

[4] C.P.Andrieux, J.M.Dumas-Bouchiat, J.M.Saveant, ‘Homogeneous Redox catalysis of electrochemical reaction,Part 

IV.Kinetic Controls in the Homogeneous Process  as characterized by Stationary and Quasi-Stationary 

Electrochemical Techniques’, J.Electroanal. Chem.,113(1980) 1-18. 

[5] C.P.Andrieux, C.Blocman, J.M.Dumas-Bouchiat, F.M.Halla, J.M.Saveant, ‘Homogeneous Redox catalysis of 

electrochemical reaction Part V. Cyclic Voltammetry’, J.Electroanal. Chem.,113(1980) 19-40. 

[6] K.S.Miller and B.Ross, ‘An Introduction to the Fractional Calculus and Fractional Differential Equations’, John 

Wiley & Sons Inc., NewYork, 1993. 

[7] Ji-Huan He, ‘Homotopy Perturbation technique’,Comput. Methods Appl. Mech. Engrg. 178(1999)    257-262. 

[8] I.Podlubny, ‘Fractional Differential Equations’, Mathematics in Science and Engineering, Volume 198 (1999). 

[9] A.C. King, J. Billingham, S.R. Otto, ‘Differential Equations: Linear, Nonlinear, Ordinary’, Partial,Cambridge 

University Press 2005. 

[10] A.A.Kilbas, H.M. Srivastava, J.J. Trujillo, ‘Theory and Applications of Fractional Differential Equations’, North-

Holland Mathematics Studies, 1st edition, Jan 12, 2006. 

[11] Ji-Huan He, ‘Homotopy Perturbation Method for solving boundary value problems’, Physics Letters A 350 (2006) 

87-88. 

[12] Dehghani.R, Ghanbari.K, Asadzadeh.M, ‘Triple Positive Solutions for Boundary value problem of a Nonlinear 

Fractional Differential Equation’, Springer, Berlin, Germany, 2007. 

[13] O.Abdulaziz, I.Hashim, S.Momani, ‘Solving system of fractional differential equations by Homotopy Perturbation 

method’, Elsevier, Physics Letters A, 372 (2008), 451-459. 

[14] Chikrii.A, Matychyn.I, ‘Riemann-Lioville, Caputo and Sequential fractional Derivative in Differential Games’, 

Birkhauser,Boston, MA, USA, 2011. 

[15] Partha S.Mallick, ‘Matlab and simulink  introduction to applications’, Fourth edition, 30 sep 2011. 

[16] Bashir.A, Nieto.J.J, ‘Sequential Fractional differential equations with three-point boundary conditions’, Comput. 

Math.Appl. 2012, 64, 3046-3052. 

[17] Ahmad El-Ajou,.O, Arqub.A, Momani.S, ‘Solving fractional two-point boundary  value problems using continuous 

analytic method’, Ain Shams Eng.J, 2013, 4, 539-547. 

[18] A.S.Yakimov, ‘Analytical solution methods for Boundary value problems’, Academic Press, 2016. 

[19] I.Ates, P.A.Zegeling, ‘A homotopy perturbation method for fractional-order advection-diffusion-reaction boundary-

value problems’, Elsevier, Applied Mathematical modelling, 47 (2017), 425-441. 

[20] N.Alghamdi, B. Ahmad, S.K. Ntouyas and A. Alsaedi, ‘Sequential fractional equations with non-local boundary 

conditions on an arbitrary interval’, Advance in Difference Equations (2017). Doi : 10.11861513662-017-1303-2. 

[21] Tariboon.J, Cuntavepanit.A, Ntouyas.S.K, Nithiarayaphaks.W, ‘Separated Boundary value problems of Sequential 

Caputo and Hadamard Fractional Differential equations’, J.Funct. Spaces 2018, 6974046. 

[22] A.S.Vatsala, G.Pageni, V.A.Vijesh, ‘Analysis of Sequential Caputo Fractional Differenial Equations versus Non-

Sequential Caputo Fractional Differential Equations with Applications’, foundations (2022), 2, 1129-1142. 

[23] Xiao-Ping Li, Et.al., ‘Application of piecewise fractional differential equation to COVID-19 infection dynamics’, 

Elsevier, Results in Physics 39(2022), 105685. 

[24] Sadia Arshad, Imran Siddique, Et.al., ‘Dynamics of a fractional order mathematical model for COVID-19 epidemic 

transmission’, Elsevier, Physica A 609(2023), 128383.