Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 5s (2024) 

 

552 
https://internationalpubls.com 

Stability and Bifurcation in a Second-Order Difference System 

 

Smitha Mary Mathew 𝟏, D.S.Dilip 𝟐, Sibi C. Babu 𝟑 

 1,2,3Department of Mathematics, St.John’s College, Anchal, Kollam District, Kerala, India. 

 1Email: smathew11@gmail.com, 2Email: dilip@stjohns.ac.in, 3Email: sibicbabu@stjohns.ac.in 

 

Article History: 

Received: 15-05-2024 

Revised: 18-06-2024 

Accepted: 11-07-2024 

Abstract:  

 In the paper, we examine the system  

 𝑥𝑛+1 = 𝛼1 + 𝑎1𝑒−𝑥𝑛−1 + 𝑏1𝑦𝑛𝑒−𝑦𝑛−1 + 𝑐1𝑒−𝑧𝑛−1 ,                       

 𝑦𝑛+1 = 𝛼2 + 𝑎2𝑒−𝑦𝑛−1 + 𝑏2𝑧𝑛𝑒−𝑧𝑛−1 + 𝑐2𝑒−𝑥𝑛−1 ,                        

 (1) 

  𝑧𝑛+1 = 𝛼3 + 𝑎3𝑒−𝑧𝑛−1 + 𝑏3𝑥𝑛𝑒−𝑥𝑛−1 + 𝑐3𝑒−𝑦𝑛−1 ,    𝑛 = 0,1,2, …, 

where 𝛼1, 𝛼2, 𝛼3, 𝑎1, 𝑎2, 𝑎3, 𝑏1, 𝑏2, 𝑏3, 𝑐1, 𝑐2, 𝑐3 are positive real numbers and the initial 

conditions 𝑥−1, 𝑥0, 𝑦−1, 𝑦0, , 𝑧−1, 𝑧0 are arbitrary nonnegative numbers. We investigate the 

persistence, boundedness, convergence, invariance, and global asymptotic character of 

the positive solutions of (1). Bifurcation diagrams are then plotted to visualize the 

periodic character. 

Keywords: persistence, boundedness, invariance, local property, global property, 

bifurcation. 

AMS Subject Classification 2000: 39A22.  

1. Introduction 

In the study of dynamical systems, difference equations play a crucial role in modeling various 

phenomena across diverse scientific disciplines, including biology, economics, engineering, and 

physics.(See [16],[21],[22],[27], [31], [32]). Unlike differential equations, which describe continuous 

change, difference equations are discrete analogues that characterize systems evolving in distinct time 

steps.(See [6],[8],[13],[22],[23],[24]) Most of the popular models like SIRS and SEIRS are mainly of 

order one. Two species models are examined in [3], [11],[29] and [33]. Competition models of two 

species second order with exponents are analyzed in [10],[12] and [15]-[20]. In [5], the authors 

analyzed a food-chain model using a first order system with three variables. More first order system 

with three variables can be seen in [1],[2],[4], [5], [7] and [28] whereas [31] deals with second order 

systems with three variables. 

The system (1) which we investigate is an extension of [30] where we focus on a system of three 

interdependent difference equations involving twelve parameters and three variables. We analyze the 

boundedness, persistence, invariance and convergence of the solutions of (1). We then plot few 

bifurcation diagrams to observe the periodic nature of the system. By exploring a three variable second 

order system, this work aims to contribute to the broader understanding of multi-parameter, multi-

variable difference system. The increased number of parameters provide a high level of flexibility to 

model real world scenarios, which is also cructial for system control dynamics. 

 

mailto:smathew11@gmail.com
mailto:dilip@stjohns.ac.in
mailto:sibicbabu@stjohns.ac.in


Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 5s (2024) 

 

553 
https://internationalpubls.com 

2. Main Results 

Theorem 2.1  The positive solution (𝑥𝑛, 𝑦𝑛, 𝑧𝑛) of (1) persists. 

It is bounded whenever  

 𝐵 = 𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3 < 1. (2) 

Proof: 

Clearly, the system persists because of the presence of 𝛼𝑖 . 

For 𝑛 = 4,5, …, (1) becomes  

 𝑥𝑛+1 ≤ 𝛼1 + 𝑎1𝑒−𝛼1 + 𝑐1𝑒−𝛼3 + 𝑏1𝑒−𝛼2[𝛼2 + 𝑎2𝑒−𝛼2 + 𝑐2𝑒−𝛼1 + 𝑏2𝑎2𝑒−𝛼3𝑧𝑛−1] 

substituting for 𝑧𝑛−1, we get  

 ≤ 𝐴1 + 𝐵𝑥𝑛−2, (3) 

 where 𝐴1 = 𝛼1 + 𝑎1𝑒−𝛼1 + 𝑐1𝑒−𝛼3 + 𝑏1𝛼2𝑒−𝛼2 + 𝑏1𝑎2𝑒−𝛼2−𝛼2 + 𝑏1𝑐2𝑒−𝛼2−𝛼1 +

𝑏1𝑏2𝛼3𝑒−𝛼2−𝛼3 + 𝑏1𝑏2𝑎3𝑒−𝛼2−𝛼3−𝛼3 + 𝑏1𝑏2𝑐3𝑒−𝛼2−𝛼2−𝛼3 

and 𝐵 = 𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3 . 

Similarly,  

 𝑦𝑛+1 ≤ 𝐴2 + 𝐵𝑦𝑛−2, (4) 

 where 𝐴2 = 𝛼2 + 𝑎2𝑒−𝛼2 + 𝑐2𝑒−𝛼1 + 𝑏2𝛼3𝑒−𝛼3 + 𝑏2𝑎3𝑒−𝛼3−𝛼3 + 𝑏2𝑐3𝑒−𝛼2−𝛼3 +

𝑏2𝑏3𝛼1𝑒−𝛼1−𝛼3 + 𝑏2𝑏3𝑎1𝑒−𝛼1−𝛼1−𝛼3 + 𝑏2𝑏3𝑐1𝑒−𝛼3−𝛼3−𝛼1. 

Also,  

 𝑧𝑛+1 ≤ 𝐴3 + 𝐵𝑧𝑛−2, (5) 

 where 𝐴3 = 𝛼3 + 𝑎3𝑒−𝛼3 + 𝑐3𝑒−𝛼2 + 𝑏3𝛼1𝑒−𝛼1 + 𝑏3𝑎1𝑒−𝛼1−𝛼1 + 𝑏3𝑐1𝑒−𝛼3−𝛼1 +

𝑏3𝑏1𝛼2𝑒−𝛼1−𝛼2 + 𝑏1𝑏3𝑎2𝑒−𝛼1−𝛼2−𝛼2 + 𝑏3𝑏1𝑐2𝑒−𝛼1−𝛼1−𝛼2. 

Now, consider the difference equations  

 𝑢𝑛+1 = 𝐴1 + 𝐵𝑢𝑛−2, 

 𝑣𝑛+1 = 𝐴2 + 𝐵𝑣𝑛−2, 

 𝑤𝑛+1 = 𝐴3 + 𝐵𝑤𝑛−2,    𝑛 = 4,5, … (6) 

 Solution (𝑢𝑛, 𝑣𝑛, 𝑤𝑛) of (6) is of the form  

 𝑢𝑛 = 𝑟1𝐵𝑛/3 + 𝑟2𝐵𝑛/3cos(
𝑛𝜋

2
) + 𝑟3𝐵𝑛/3sin(

𝑛𝜋

2
) +

𝐴1

1−𝐵
,    𝑛 = 5,6, …, (7) 

  𝑣𝑛 = 𝑠1𝐵𝑛/3 + 𝑠2𝐵𝑛/3cos(
𝑛𝜋

2
) + 𝑠3𝐵𝑛/3sin(

𝑛𝜋

2
) +

𝐴2

1−𝐵
,    𝑛 = 5,6 …, (8) 

  𝑤𝑛 = 𝑝1𝐵𝑛/3 + 𝑝2𝐵𝑛/3cos(
𝑛𝜋

2
) + 𝑝3𝐵𝑛/3sin(

𝑛𝜋

2
) +

𝐴3

1−𝐵
,    𝑛 = 5,6 …, (9) 

 where 𝑝𝑖, 𝑟𝑖, 𝑠𝑖, 𝑖 = 1,2,3 depend on 𝑤4, 𝑢4, 𝑣4 respectively. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 5s (2024) 

 

554 
https://internationalpubls.com 

Hence, (𝑢𝑛, 𝑣𝑛, 𝑤𝑛) is bounded.  

Now we examine (𝑢𝑛, 𝑣𝑛, 𝑤𝑛) such that the initial conditions of 6 and 1 are same. Clearly we can 

conclude that (𝑥𝑛, 𝑦𝑛, 𝑧𝑛) is bounded.  

Theorem 2.2  Let (2) hold. Let 𝐴1, 𝐴2, 𝐴3 be defined as in Theorem 2.1. Then [𝛼1,
𝐴1

1−𝐵
] 

× [𝛼2,
𝐴2

1−𝐵
] × [𝛼3,

𝐴3

1−𝐵
] is an invariant set for the system (1).  

Proof: Take 𝐼1 = [𝛼1,
𝐴1

1−𝐵
], 𝐼2 = [𝛼2,

𝐴2

1−𝐵
] and 𝐼3 = [𝛼3,

𝐴3

1−𝐵
]. 

Let 𝑥−1, 𝑥0 ∈ 𝐼1, 𝑦−1, 𝑦0 ∈ 𝐼2 and 𝑧−1, 𝑧0 ∈ 𝐼3. 

Then 𝑥1 ≤ 𝛼1 + 𝑎1𝑒−𝛼1 + 𝑐1𝑒−𝛼3 + 𝑏1𝑒−𝛼2𝑦0 

Since 𝑦0 ≤
𝐴2

1−𝐵
, we get  

𝑥1 ≤ [
𝛼1 + 𝑎1𝑒−𝛼1 + 𝑐1𝑒−𝛼3 − 𝑏1𝑏2𝑏3𝛼1𝑒−𝛼1−𝛼2−𝛼3 − 𝑏1𝑏2𝑏3𝑎1𝑒−𝛼1−𝛼1−𝛼2−𝛼3

1 − 𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3
] 

+[
−𝑏1𝑏2𝑏3𝑐1𝑒−𝛼1−𝛼2−𝛼3−𝛼3 + 𝑏1𝛼2𝑒−𝛼2 + 𝑏1𝑎2𝑒−𝛼2−𝛼2 + 𝑏1𝑐2𝑒−𝛼1−𝛼2

1 − 𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3
] 

+[
𝑏1𝑏2𝛼3𝑒−𝛼2−𝛼3 + 𝑏1𝑏2𝑎3𝑒𝛼3−𝛼3 + 𝑏1𝑏2𝑐3𝑒−𝛼2−𝛼2−𝛼3 + 𝑏1𝑏2𝑏3𝛼1𝑒−𝛼2−𝛼2−𝛼3

1 − 𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3
] 

+[
𝑏1𝑏2𝑏3𝑎1𝑒−𝛼1−−𝛼1−𝛼2−𝛼3 + 𝑏1𝑏2𝑏3𝑐1𝑒−𝛼1−𝛼2−𝛼3−𝛼3

1 − 𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3
] 

i.e., 𝑥1 ∈ 𝐼1. Similarly we get 𝑦1 ∈ 𝐼2 and 𝑧1 ∈ 𝐼3. 

Hence the induction the proof follows.  

Theorem 2.3  Assume (2). Let 𝐴1, 𝐴2, 𝐴3 be as in the Theorem 2.1. Let 𝐼4 = [𝛼1,
𝐴1+𝜖

1−𝐵
], 𝐼5 = [𝛼2,

𝐴2+𝜖

1−𝐵
] 

and 𝐼6 = [𝛼3,
𝐴3+𝜖

1−𝐵
] where 𝜖 is arbitrary. Then 𝑥𝑛 ∈ 𝐼4, 𝑦𝑛 ∈ 𝐼5 and 𝑧𝑛 ∈ 𝐼6, for every 𝑛 ≥ 𝑁, 𝑁 ∈ ℕ.  

Proof: 

Given (𝑥𝑛, 𝑦𝑛, 𝑧𝑛) be a positive solution of (1). 

Theorem 2.1 implies, 

limsup𝑛→∞𝑥𝑛 = 𝑃 < ∞, limsup𝑛→∞𝑦𝑛 = 𝑄 < ∞ and limsup𝑛→∞𝑧𝑛 = 𝑅 < ∞. 

Theorem 2.1 implies, 

𝑥𝑛+1 ≤ 𝐴1 + 𝑏1𝑏2𝑏3𝑥𝑛−2𝑒−𝛼1−𝛼2−𝛼3 , 𝑦𝑛+1 ≤ 𝐴2 + 𝑏1𝑏2𝑏3𝑦𝑛−2𝑒−𝛼1−𝛼2−𝛼3 

and 

𝑧𝑛+1 ≤ 𝐴3 + 𝑏1𝑏2𝑏3𝑧𝑛−2𝑒−𝛼1−𝛼2−𝛼3 . 

Hence, 𝑃 ≤
𝐴1

1−𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3
, 𝑄 ≤

𝐴2

1−𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3
  and 𝑅 ≤

𝐴3

1−𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3
. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 5s (2024) 

 

555 
https://internationalpubls.com 

Hence, the result.  

Here we state a lemma which is an an extension of Lemma 5 in [10] and a variation of Theorem 1.16 

in [26].  

Lemma 2.4  Assume 𝐴, 𝐵, 𝐶, 𝐷, 𝐸, 𝐹 represent reals. Let 𝑓1: [𝐴, 𝐵] × [𝐶, 𝐷] × [𝐶, 𝐷] × [𝐸, 𝐹] →

[𝐴, 𝐵], 𝑓2: [𝐶, 𝐷] × [𝐸, 𝐹] × [𝐸, 𝐹] × [𝐴, 𝐵] → [𝐶, 𝐷] and 𝑓3: [𝐸, 𝐹] × [𝐴, 𝐵] × [𝐴, 𝐵] × [𝐶, 𝐷] →

[𝐸, 𝐹] be continuous. Examine  

 𝑥𝑛+1 = 𝑓1(𝑥𝑛−1, 𝑦𝑛, 𝑦𝑛−1, 𝑧𝑛−1),                 

 𝑦𝑛+1 = 𝑓2(𝑦𝑛−1, 𝑧𝑛, 𝑧𝑛−1, 𝑥𝑛−1),                 (10) 

 𝑧𝑛+1 = 𝑓3(𝑧𝑛−1, 𝑥𝑛, 𝑥𝑛−1, 𝑦𝑛−1),    𝑛 = 0,1,2, … 

where 𝑥−1, 𝑥0 ∈ [𝐴, 𝐵], 𝑦−1, 𝑦0 ∈ [𝐶, 𝐷] and 𝑧−1, 𝑧0 ∈ [𝐸, 𝐹]. (or 𝑥𝑛0
, 𝑥𝑛0+1 ∈ [𝐴, 𝐵], 𝑦𝑛0

, 𝑦𝑛0+1 ∈

[𝐶, 𝐷], 𝑧𝑛0
, 𝑧𝑛0+1 ∈ [𝐸, 𝐹], 𝑛0 ∈ ℕ). Assume the conditions given below holds.   

    1.  If 𝑓1(𝑥, 𝑦, 𝑧, 𝑢), 𝑓2(𝑥, 𝑦, 𝑧, 𝑢) and 𝑓3(𝑥, 𝑦, 𝑧, 𝑢) are nonincreasing in x, nondecreasing in y, 

nonincreasing in z and nonincreasing in u.  

    2.  If (𝑚1, 𝑀1, 𝑚2, 𝑀2, 𝑚3, 𝑀3) ∈ [𝐴, 𝐵]2 × [𝐶, 𝐷]2 × [𝐸, 𝐹]2 satisfies the systems 𝑚1 =

𝑓1(𝑀1, 𝑚2, 𝑀2, 𝑀3); 𝑀1 = 𝑓1(𝑚1, 𝑀2, 𝑚2, 𝑚3), 𝑚2 = 𝑓2(𝑀2, 𝑚3, 𝑀3, 𝑀1); 𝑀2 = 𝑓2(𝑚2, 𝑀3, 𝑚3, 𝑚1) 

and 𝑚3 = 𝑓3(𝑚1, 𝑀1, 𝑀3, 𝑀2); 𝑀3 = 𝑓3(𝑀1, 𝑚1, 𝑚3, 𝑚2) then 𝑚1 = 𝑀1, 𝑚2 = 𝑀2 and 𝑚3 = 𝑀3,  

then (�̅�, �̅�, 𝑧̅) is the unique equilibrium point of (10) where �̅� ∈ [𝐴, 𝐵],�̅� ∈ [𝐶, 𝐷] and 𝑧̅ ∈ [𝐸, 𝐹]. And 

any other solution of (10) converges to (�̅�, �̅�, 𝑧̅).  

Theorem 2.5  Let (2) hold. Suppose  

 𝑎1𝑒−𝛼1 < 1, 𝑎2𝑒−𝛼2 < 1, 𝑎3𝑒−𝛼3 < 1 (11) 

 and  

 
[𝐷2𝐷3+𝐵3𝐿2][𝐷1𝐷2+𝐵2𝐿1][𝐷3𝐷1+𝐵1𝐿3]

[𝐵2𝐵3−𝐷2𝐿3][𝐵1𝐵2−𝐷1𝐿2][𝐵3𝐵1−𝐷3𝐿1]
< 1, (12) 

 where 𝐵1 = 1 − 𝑎1𝑒−𝛼1 , 𝐵2 = 1 − 𝑎2𝑒−𝛼2 , 𝐵3 = 1 − 𝑎3𝑒−𝛼3 , 𝐷1 = 𝑏1𝑒−𝛼2(1 +
𝐴2

1−𝐵
), 𝐷2 =

𝑏2𝑒−𝛼3(1 +
𝐴3

1−𝐵
), 𝐷3 = 𝑏3𝑒−𝛼1(1 +

𝐴1

1−𝐵
), 𝐿1 = 𝑐1𝑒−𝛼3 , 𝐿2 = 𝑐2𝑒−𝛼1 , 𝐿3 = 𝑐3𝑒−𝛼2 . Then 𝐸(�̅�, �̅�, 𝑧̅) is 

the unique positive equilibrium of (1). And any solution of (1) converges to 𝐸(�̅�, �̅�, 𝑧̅).  

Proof: 

Define                  𝑓1(𝑥, 𝑦, 𝑧) = 𝛼1 + 𝑎1𝑒−𝑥 + 𝑏1𝑦𝑒−𝑦 + 𝑐1𝑒−𝑧,  

𝑓2(𝑥, 𝑦, 𝑧) = 𝛼2 + 𝑎2𝑒−𝑦 + 𝑏2𝑧𝑒−𝑧 + 𝑐2𝑒−𝑥,        𝑓3(𝑥, 𝑦, 𝑧) = 𝛼3 + 𝑎3𝑒−𝑧 + 𝑏3𝑥𝑒−𝑥 + 𝑐3𝑒−𝑦𝑆. 

Take 𝑚𝑖 ≤ 𝑀𝑖, 𝑖 = 1,2,3 to denote positive reals where and  

𝑚1 = 𝛼1 + 𝑎1𝑒−𝑀1 + 𝑏1𝑚2𝑒−𝑀2 + 𝑐1𝑒−𝑀3 , 𝑀1 = 𝛼1 + 𝑎1𝑒−𝑚1 + 𝑏1𝑀2𝑒−𝑚2 + 𝑐1𝑒−𝑚3 , 

 

 𝑚2 = 𝛼2 + 𝑎2𝑒−𝑀2 + 𝑏2𝑚3𝑒−𝑀3 + 𝑐1𝑒−𝑀1 , 𝑀2 = 𝛼2 + 𝑎2𝑒−𝑚2 + 𝑏2𝑀3𝑒−𝑚3 + 𝑐1𝑒−𝑚1 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 5s (2024) 

 

556 
https://internationalpubls.com 

and  

 𝑚3 = 𝛼3 + 𝑎3𝑒−𝑀3 + 𝑏3𝑚1𝑒−𝑀3 + 𝑐1𝑒−𝑀2 , 𝑀3 = 𝛼3 + 𝑎3𝑒−𝑚3 + 𝑏3𝑀1𝑒−𝑚1 + 𝑐1𝑒−𝑚2 . (13) 

 Therefore, 𝑀1 − 𝑚1 = 𝑎1[𝑒−𝑚1 − 𝑒−𝑀1] + 𝑏1[𝑀2𝑒−𝑚2 − 𝑚2𝑒−𝑀2] + 𝑐1[𝑒−𝑚3 − 𝑒−𝑀3].  

 𝑀1 − 𝑚1 = 𝑎1[𝑒−𝑚1 − 𝑒−𝑀1] + 𝑏1𝑒−𝑚2−𝑀2[𝑀2𝑒𝑀2 − 𝑚2𝑒𝑚2] + 𝑐1[𝑒−𝑚3 − 𝑒−𝑀3]. (14) 

 Here, there exists a 𝜁1 , 𝑀2 ≥ 𝜁1 ≥ 𝑚2 satisfying  

 𝑀2𝑒𝑀2 − 𝑚2𝑒𝑚2 = (1 + 𝜁1)𝑒1
𝜁

(𝑀2 − 𝑚2). (15) 

 From (14) and (15) we get,  

 𝑀1 − 𝑚1 = 𝑎1[𝑒−𝑚1 − 𝑒−𝑀1] + 𝑏1𝑒−𝑚2−𝑀2+𝜁1(1 + 𝜁1)[𝑀2 − 𝑚2] + 𝑐1[𝑒−𝑚3 − 𝑒−𝑀3]. (16) 

 Now, 𝑎1[𝑒−𝑚1 − 𝑒−𝑀1] = 𝑎1𝑒−𝑚1−𝑀1[𝑒𝑀1 − 𝑒𝑚1]. 

And, there exists a 𝜆, 𝑀1 ≥ 𝜆 ≥ 𝑚1 satisfying  

 𝑎1[𝑒−𝑚1 − 𝑒−𝑀1] = 𝑎1𝑒−𝑚1−𝑀1+𝜆[𝑀1 − 𝑚1]. (17) 

 Since 𝑚1, 𝑀1 ≥ 𝛼1,  

 𝑎1[𝑒−𝑚1 − 𝑒−𝑀1] ≤ 𝑎1𝑒−𝛼1[𝑀1 − 𝑚1]. (18) 

 Thus from (16) and (18) we get,  

 𝑀1 − 𝑚1 ≤ 𝑎1𝑒−𝛼1[𝑀1 − 𝑚1] + 𝑏1𝑒−𝑚2−𝑀2+𝜁1(1 + 𝜁1)[𝑀2 − 𝑚2] + 𝑐1𝑒−𝛼3[𝑀3 − 𝑚3]. (19) 

 Since 𝑚2, 𝑀2 ≥ 𝛼2, (19) becomes  

 𝑀1 − 𝑚1 ≤ 𝑎1𝑒−𝛼1[𝑀1 − 𝑚1] + 𝑏1𝑒−𝛼2(1 + 𝜁1)[𝑀2 − 𝑚2] + 𝑐1𝑒−𝛼3[𝑀3 − 𝑚3]. (20) 

 i.e.,  

 [1 − 𝑎1𝑒−𝛼1][𝑀1 − 𝑚1] ≤ 𝑏1𝑒−𝛼2(1 + 𝜁1)[𝑀2 − 𝑚2] + 𝑐1𝑒−𝛼3[𝑀3 − 𝑚3]. (21) 

Also, (13) can be written as  

 𝑀2 = 𝛼2 + 𝑎2𝑒−𝑚2 + 𝑏2[𝛼3 + 𝑎3𝑒−𝑚3 + 𝑏3𝑀1𝑒−𝑚1 + 𝑐3𝑒−𝑚2]𝑒−𝑚3 + 𝑐2𝑒−𝑚1 . (22) 

Substituting again for 𝑀1 and simplifying we get  

 𝑀2 ≤
𝐴2

1−𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3
. (23) 

 Since 𝜁1 ≤ 𝑀2 we get,  

 𝜁1 ≤
𝐴2

1−𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3
. (24) 

 Therefore, (21) becomes  

 [1 − 𝑎1𝑒−𝛼1][𝑀1 − 𝑚1] ≤ 𝑏1𝑒−𝛼2[1 +
𝐴2

1−𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3
][𝑀2 − 𝑚2] + 𝑐1𝑒−𝛼3[𝑀3 − 𝑚3].(25) 

 Similarly we get,  

 [1 − 𝑎2𝑒−𝛼2][𝑀2 − 𝑚2] ≤ 𝑏2𝑒−𝛼3[1 +
𝐴3

1−𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3
][𝑀3 − 𝑚3] + 𝑐2𝑒−𝛼1[𝑀1 − 𝑚1](26) 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 5s (2024) 

 

557 
https://internationalpubls.com 

 and  

 [1 − 𝑎3𝑒−𝛼3][𝑀3 − 𝑚3] ≤ 𝑏3𝑒−𝛼1[1 +
𝐴1

1−𝑏1𝑏2𝑏3𝑒−𝛼1−𝛼2−𝛼3
][𝑀1 − 𝑚1] + 𝑐3𝑒−𝛼2[𝑀2 − 𝑚2](27) 

 From (25), (26) and (27) we get, 

[𝐵1𝐵2 − 𝐷1𝐿2]

𝐵2
[𝑀1 − 𝑚1] ≤

[𝐷1𝐷2 + 𝐵2𝐿1]

𝐵2
[𝑀3 − 𝑚3]. 

Similarly,  

 
[𝐵2𝐵3−𝐷2𝐿3]

𝐵3
[𝑀2 − 𝑚2] ≤

[𝐷2𝐷3+𝐵3𝐿2]

𝐵3
[𝑀1 − 𝑚1] (28) 

 and  

 
[𝐵3𝐵1−𝐷3𝐿1]

𝐵1
[𝑀3 − 𝑚3] ≤

[𝐷3𝐷1+𝐵1𝐿3]

𝐵1
[𝑀2 − 𝑚2] (29) 

 Hence from (12) and (28), we get 𝑀1 = 𝑚1. 

Similarly we get 𝑀2 = 𝑚2 and 𝑀3 = 𝑚3.  

Hence by Lemma 2.4, we get the required result.  

 

Theorem 2.6  Assume (2), (11) and (12) hold. If 

𝑎1𝑒−𝛼1[1 + 𝑎2𝑒−𝛼2] + 𝑎2𝑒−𝛼2[1 + 𝑎3𝑒−𝛼3] + 𝑎3𝑒−𝛼3[1 + 𝑎1𝑒−𝛼1] + 𝑎1𝑎2𝑎3𝑒−𝛼1𝑒−𝛼2𝑒−𝛼3 

+
𝐵

(1−𝐵)3 [(1 − 𝐵)3 + (1 − 𝐵)2(𝐴1 + 𝐴2 + 𝐴3)  

 +(1 − 𝐵)(𝐴1𝐴2 + 𝐴2𝐴3 + 𝐴1𝐴3) + 𝐴1𝐴2𝐴3] < 1, (30) 

 where 𝐵, 𝐴1, 𝐴2, 𝐴3 are as in Theorem 2.1, then 𝐸(�̅�, �̅�, 𝑧̅) of (2.5) is globally asymptotically stable.  

 Proof: 

We need to illustrate 𝐸(�̅�, �̅�, 𝑧̅) is locally asymptotic. Construct the Jacobian 𝐽𝐹(�̅�, �̅�, 𝑧̅) about 

𝐸(�̅�, �̅�, 𝑧̅). Its characteristic equation is given by 

𝜆6 + 𝜆4(𝑎1𝑒−�̅� + 𝑎2𝑒−�̅� + 𝑎3𝑒−�̅�) + 𝜆3(𝑏2𝑐3𝑒−�̅�−�̅� + 𝑏1𝑐2𝑒−�̅�−�̅� + 𝑏3𝑐1𝑒−�̅�−�̅� + 𝑏1𝑏2𝑏3𝑒−�̅�−�̅�−�̅�)

+ 𝜆2(𝑎2𝑎3𝑒−�̅�−�̅� − 𝑏2𝑐3𝑧̅𝑒−�̅�−�̅� + 𝑎1𝑎3𝑒−�̅�−�̅� − 𝑏3𝑐1�̅�𝑒−�̅�−�̅� + 𝑎1𝑎2𝑒−�̅�𝑒−�̅�

− 𝑏1𝑐2�̅�𝑒−�̅�𝑒−�̅� + 𝑒−�̅�−�̅�−�̅�[𝑏1𝑏2𝑏3�̅� + 𝑏1𝑏2𝑏3�̅� + 𝑏1𝑏2𝑏3𝑧̅])

+ 𝜆𝑒−�̅�−�̅�−�̅�(−𝑏1𝑏2𝑏3�̅��̅� − 𝑏1𝑏2𝑏3�̅�𝑧̅ − 𝑏1𝑏2𝑏3�̅�𝑧̅ + 𝑎3𝑏1𝑐2 + 𝑎1𝑏2𝑐3 + 𝑎2𝑏3𝑐1)

+ 𝑒−�̅�−�̅�−�̅�(𝑎1𝑎2𝑎3 + 𝑏1𝑏2𝑏3�̅��̅�𝑧̅ + 𝑐1𝑐2𝑐3 − 𝑎2𝑏3𝑐1�̅� − 𝑎3𝑏1𝑐2�̅� − 𝑎1𝑏2𝑐3𝑧̅) = 0. 

Applying Remark 1.3.1 of [25], 

|𝑎1𝑒−�̅�| + |𝑎2𝑒−�̅�| + |𝑎3𝑒−�̅�| + |𝑏1𝑏2𝑏3𝑒−�̅�−�̅�−�̅�| + |𝑏1𝑏2𝑏3�̅�𝑒−�̅�−�̅�−�̅�| + |𝑏1𝑏2𝑏3�̅�𝑒−�̅�−�̅�−�̅�|

+ |𝑏1𝑏2𝑏3𝑧̅𝑒−�̅�−�̅�−�̅�| + |𝑎1𝑎2𝑒−�̅�𝑒−�̅�| + |𝑎1𝑎3𝑒−�̅�𝑒−�̅�| + |𝑎2𝑎3𝑒−�̅�𝑒−�̅�|

+ |𝑏1𝑏2𝑏3�̅��̅�𝑒−�̅�−�̅�−�̅�| + |𝑏1𝑏2𝑏3�̅�𝑧̅𝑒−�̅�−�̅�−�̅�| + |𝑏1𝑏2𝑏3�̅�𝑧̅𝑒−�̅�−�̅�−�̅�|

+ |𝑎1𝑎2𝑎3𝑒−�̅�−�̅�−�̅�| + |𝑏1𝑏2𝑏3�̅��̅�𝑧̅𝑒−�̅�−�̅�−�̅�| < 1 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 5s (2024) 

 

558 
https://internationalpubls.com 

is satisfied whenever 𝑎1𝑒−𝛼1[1 + 𝑎2𝑒−𝛼2] + 𝑎2𝑒−𝛼2[1 + 𝑎3𝑒−𝛼3] + 𝑎3𝑒−𝛼3[1 + 𝑎1𝑒−𝛼1] +  

 𝑎1𝑎2𝑎3𝑒−𝛼1𝑒−𝛼2𝑒−𝛼3 + 𝐵[�̅��̅�𝑧̅ + �̅��̅� + �̅�𝑧̅ + �̅�𝑧̅ + 𝑧̅ + �̅� + �̅� + 1] < 1. (31) 

 Clearly from Theorem 2.1, 

 �̅� ≤
𝐴1

1−𝐵
, (32) 

 �̅� ≤
𝐴2

1−𝐵
 (33) 

 and  

 𝑧̅ ≤
𝐴3

1−𝐵
. (34) 

 Substitute (32), (33) and (34) in (31). 

Use Remark 1.3.1 of [25] and Theorem 2.5 to get the result.  

3. Numerical Analysis and Open Problem 

In this section we observe the dynamics of the discrete model (1) numerically and propose an open 

problem. Figure (1a) shows the bifurcation diagram with 𝑎1 as bifurcation parameter and figure (1b) 

shows the plots of 𝑥𝑛, 𝑦𝑛, 𝑧𝑛 for a particular value,ie., 𝑎1 = 32.0. Figure (1b) shows that the plot is 

eventually 4-periodic. Similarly figures(2a) and (3a) shows the bifurcation diagrams with 𝑎2 and 𝑎3 

as bifurcation parameter, whereas figures(2b) and (3b) shows their corresponding plots for 𝑎2 = 10.0 

and 𝑎3 = 45.9 respectively. Figures (2b) and (3b) shows that the plots are eventually 4-periodic.  

3.1 Open Problem 

Derive the condition for (1) to be eventually 4-periodic. 

 

a) [Bifurcation Diagrams of (1) with 𝑎1as  bifurcation parameter]                    b)  [Plots of 𝑥𝑛, 𝑦𝑛, 𝑧𝑛 

with 𝑎1 = 32.0] 

Figure  1: Here 𝛼1 = 2.3, 𝛼2 = 3.2, 𝛼3 = 2.6, 𝑏1 = .4, 𝑐1 = .5, 𝑎2 = 0.5, 𝑏2 = 4.4, 𝑐2 = 3.5, 𝑎3 =

0.9, 𝑏3 = 0.6, 𝑐3 = 0.4. 

 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 5s (2024) 

 

559 
https://internationalpubls.com 

 

a) [Bifurcation Diagrams of (1) with 𝑎2 as  

bifurcation parameter] 

b)  [Plots of 𝑥𝑛, 𝑦𝑛, 𝑧𝑛 with 𝑎2 = 10.0] 

Figure  2: Here 𝛼1 = 5.3, 𝛼2 = 0.2, 𝛼3 = 12.6, 𝑎1 = 0.5, 𝑏1 = .4, 𝑐1 = .5, 𝑏2 = 4.4, 𝑐2 = 1.5, 𝑎3 =

5.9, 𝑏3 = 8.6, 𝑐3 = 0.4. 

 

a) [Bifurcation Diagrams of (1) with 𝑎3 as  

bifurcation parameter] 

b)  [Plots of 𝑥𝑛, 𝑦𝑛, 𝑧𝑛 with 𝑎3 = 45.0] 

Figure  3: Here 𝛼1 = 2.3, 𝛼2 = 3.2, 𝛼3 = 2.6, 𝑎1 = 0.9, 𝑏1 = .4, 𝑐1 = .5, 𝑎2 = 0.5, 𝑏2 = 4.4, 𝑐2 =

3.5, 𝑏3 = 0.6, 𝑐3 = 0.4. 

4. Conclusion 

In this paper, we studied the dynamics of a second-order system defined by three variables, focusing 

on the existence of a unique positive equilibrium and its global stability. This study is particularly 

relevant in the context of population biology, where understanding the conditions for local asymptotic 

stability and global stability are crucial. We successfully established the conditions which assure the 

global asymptotic stability of the unique positive equilibrium.Moreover, we proposed an open problem 

that invite further investigation into the conditions necessary for the system to exhibit 4-periodic 

behavior. Addressing this problem will provide deeper insights into the periodic nature of the system 

and its long-term behavior.  



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 5s (2024) 

 

560 
https://internationalpubls.com 

References 

[1] A. S. Ackleh and P. Zhang, Competitive Exclusion in a Discrete Stage-Structured Two Species Model, Mathematical 

Modelling of Natural Phenomena, 4(6), (2009), 156 - 175.  

[2] Azmy S. Ackleh and Patrick De Leenheer, Discrete three-stage population model: persistence and global stability 

results, Journal of Biological Dynamics, 2(4), (2008), 415 - 427. 

[3] A. T. Ademola, P. O. Arawomo, and A. S. Idowu,Stability, boundedness and periodic solutions to certain second 

order delay differential equations, Proyecciones (Antofagasta, On line), vol. 36, no. 2, pp. 257-282, Jun. 2017.  

[4] Irfan Ali, Umer Saeed, Qamar Din, Bifurcation analysis and chaos control in discrete-time system of three competing 

species, Arabian Journal of Mathematics, 8, (2019), 1-14.  

[5] Ll. Alseda, B. Vidiella, R. Sole, J.T. Lazaro, J. Sardanyes, Dynamics in a ́ time-discrete food-chain model with strong 

pressure on preys, Communications in Nonlinear Science and Numerical Simulation, 84, (2020), 105187. 

[6] J. Leo Amalraj, M.Maria Susai Manuel, Dumitru Baleanu and D.S.Dilip, Global Stability, Periodicity and Bifurcation 

Analysis of a Difference Equation, AIP Advances, 13, (2023), 015116.  

[7] Ritwick Banerjee, Pritha Das, Debasis Mukherjee, Stability and permanence of a discrete-time two-prey one-predator 

system with Holling Type-III functional response, Chaos, Solitons and Fractals, 117, (2001), 240 - 248.  

[8] C.A. Clark, M. R. S Kulenovic and James F. Selgrade, On a Sytem of Rational Difference Equations, Journal of 

Difference Equations and Applications, 11(7), (2005), 565 - 580. 

[9] Elias Camouzis and G. Ladas, Dynamics of Third Order Rational Difference Equations with Open Problems, 

Chapman & Hall/CRC, 2007. 

[10] D S Dilip and Smitha Mary Mathew, Dynamics of a Second Order Nonlinear Difference System with Exponents, 

Journal of the Egyptian Mathematical Society, (2021) 29:10. 

[11] D S Dilip and Smitha Mary Mathew, Stability analysis of a time varying population model without migration, Journal 

of Difference Equations and Applications, 27(11), (2021), 1525-1536. 

[12] D.S.Dilip and Tony Philip, Stability and Bifurcation Analysis of Two Spatial Population Dynamics Models, 

International Journal of Applied and Computational Mathematics, 8, 108(2022). 

[13] N. Fotiades and G. Papaschinopoulos, Existence, Uniqueness and Attractivity of Prime Period Two Solution for a 

Difference Equation of Exponential form, Applied Mathematics and Computation, 218 (2012), 11648 - 11653. 

[14] M. Kulenovic and G. Ladas, Dynamics of Second Order Rational Difference Equations, Chapman & Hall/CRC, 2002.  

[15] G. Papaschinopoulos, M. A. Radin and C.J. Schinas, On the system of Two Difference Equations of Exponential Form, 

Mathematical and Computer Modelling, 54 (2011), 2969 - 2977.  

[16] G. Papaschinopoulos and C.J. Schinas, On the Dynamics of two Exponential Type Systems of Difference Equations, 

Computers and Mathematics with Applications, 64 (2012), 2326 - 2334.  

[17] G. Papaschinopoulos, G. Ellina and K.B.Papadopoulos, Asymptotic Behavior of the Positive Solutions of an 

Exponential Type System of Difference Equations, Applied Mathematics and Computation, 245 (2014), 181 - 190. 

[18] G. Papaschinopoulos, C.J. Schinas and G. Ellina, On the Dynamics of the Solutions of a Biological Model, Journal of 

Difference Equations and Applications, 20(5 - 6), (2014), 694 - 705.  

[19] G. Papaschinopoulos, N. Fotiades and C.J. Schinas, On a System of Difference Equations Including Negative 

Exponential Terms, Journal of Difference Equations and Applications, 20(5 - 6), (2014), 717 - 732.  

[20] N. Psarros and G. Papaschinopoulos, Long-term Behavior of Positive Solutions of an Exponentially Self-regulating 

System of Difference Equations, International Journal of Biomathematics, 10(3), (2017), 1750045. 

[21] Muhammad Naeem Qureshi, A. Qadeer Khan and Qamar Din, Asymptotic Behavior of a Nicholson-Bailey Model, 

Advances in Difference Equations, (2014), 2014:62. 

[22] Samra Moranjkic and Zehra Nurkanovic, Basins of attraction of certain rational anti-competitive system of difference 

equations in the plane, Advances in Difference Equations, (2012), 2012:153. 

[23] Hui Feng, Huili Ma and Wandi Ding, Global Asymptotic Behavior of Positive Solutions for Exponential Form 

Difference Equations with Three Parameters, Journal of Applied Analysis and Computation, 6(3), (2016), 600 - 606. 

[24] Q Din, MN Qureshi and A Qadeer Khan, Dynamics of a fourth-order system of rational difference equations, 

Advances in Difference Equations, (2012), 2012:215. 

[25] Kocic VL and Ladas G, Global Behavior of Nonlinear Difference Equations of Higher Order with Applications, 

Kluwer Academic Publishers: Dordrecht/Boston/London, 1993. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 31 No. 5s (2024) 

 

561 
https://internationalpubls.com 

[26] E.A.Grove and Ladas G, Periodicities in Nonlinear Difference Equations, Chapman & Hall/CRC, 2005. 

[27] Qamar Din, Complexity and Chaos Control in a Discrete-Time Prey-Predator Model, Communications in Nonlinear 

Science and Numerical Simulation, 49 (2017), 113 - 134.  

[28] Guichen Lu and Zhengyi Lu, Non-permanence for three-species Lotka-Volterra cooperative difference systems, 

Advances in Difference Equations, (2017), 2017:152. 

[29] Sibi C. Babu, D.S.Dilip and Smitha Mary Mathew, Behavior of solutions of a discrete population model with 

mutualistic interaction, Computational and Mathematical Biophysics, 12(1), (2024), 20230121. 

[30] Smitha Mary Mathew and D.S.Dilip, Dynamics of a second order three species nonlinear difference system with 

exponents, Proyecciones (Antofagasta, On line), 41(4), (2022), 983-997. 

[31] Smitha Mary Mathew and D.S.Dilip, Dynamics of interspecific 𝑘 species competition model, Journal of 

Interdisciplinary Mathematics, 25(3), (2022), 629-638.  

[32] D. Tilman and D. Wedin, Oscillations and Chaos in the Dynamics of a Perennial Grass, Letters to Nature, 353 (1991), 

653 - 655. 

[33] Yu X., Zhu Z. and Li Z, Stability and bifurcation analysis of two-species competitive model with Michaelis–Menten 

type harvesting in the first species, Adv Differ Equ 2020, 397 (2020). https://doi.org/10.1186/s13662-020-02817-4.