Int. J. Anal. Appl. (2022), 20:38 Exact Solutions and Stability of Fourth Order Systems of Difference Equations Using Padovan Numbers Marwa M. Alzubaidi∗ Department of Mathematics, College of Duba, University of Tabuk, Saudi Arabia ∗Corresponding author: mmialzubaidi@hotmail.com Abstract. Difference equations are widely utilized to describe some phenomena arising in nonlinear sciences. In particular, systems of difference equations play an important role in investigating most nonlinear applications. Future behaviors of such phenomena can be sometimes known and understood by using exact solutions of systems of difference equations. Therefore, this article investigates the exact solutions of fourth order systems of difference equations. We use successive iterations and Padovan numbers to obtain the exact solutions in the form of rational functions. The stability of the considered systems are analyzed using Jacobian matrix. Real equilibrium points are found saddle. Under some selected parameters, we plot some 2D figures to show the behavior of the obtained solutions. The used methods can be successfully applied for high order systems of difference equations. 1. Introduction Differential equations have been widely used to model various economical, physical, biological and artificial phenomena. These equations describe populations or objects in which time is continuous. In contrast, difference equations have been extensively used to describe populations or objects that evolve discrete time. The study of the theory of difference equations has strongly become an active topic for some researchers. The main reason behind that is that difference equations are used to model and describe most natural and non-natural phenomena such as those occurred in physics, biology, chemistry, economy, engineering, etc. Difference equations appear as discrete mathematical models of these phenomena. Recursive equations model wide spectrum of biomedical phenomena such as cell proliferation, cancer growth and genetics [1]. Discrete dynamical systems are also used to describe Received: Jul. 3, 2022. 2010 Mathematics Subject Classification. 39A30. Key words and phrases. equilibrium; stability; system of difference equations; Padovan numbers; plastic number. https://doi.org/10.28924/2291-8639-20-2022-38 ISSN: 2291-8639 © 2022 the author(s). https://doi.org/10.28924/2291-8639-20-2022-38 2 Int. J. Anal. Appl. (2022), 20:38 various applications in different sciences. Moreover, difference equations play a key role in investigating the solutions of various differential equations. The solutions and the stability of many systems of difference equations have been recently discussed. For instance, Tollu et al. [2] used Fibonacci numbers to investigate the solutions of the difference equations xn+1 = 1 1 + xn , yn+1 = 1 −1 + yn . Yazlik et al. [3] obtained the solutions of the following systems xn+1 = xn−1 ± 1 ynxn−1 , yn+1 = yn−1 ± 1 xnyn−1 . (1.1) The author in [4] investigated the solutions of the the system xn+1 = xn−1yn−3 yn−1 (−1 −xn−1yn−3) , yn+1 = yn−1xn−3 xn−1 (±1 ±yn−1xn−3) , n = 0, 1, ..., Moreover, Almatrafi and Alzubaidi [5] investigated the stability, periodicity and the solutions of the difference equation xn+1 = c1xn−3 + c2xn−3 c3xn−3 −c4xn−7 . Cinar [6] obtained the positive solutions of the following system xn+1 = 1 yn , yn+1 = yn yn−1xn−1 . In [7], Elsayed extracted the solutions of the following system of difference equations xn+1 = xn−1 ±1 + xn−1yn , yn+1 = yn−1 ∓1 + yn−1xn . Alayachi et al. [8] discussed the solutions of nonlinear rational systems of difference equations of third order given by xn+1 = ynzn−1 yn ±xn−2 , yn+1 = znxn−1 zn ±yn−2 , zn+1 = xnyn−1 xn ±zn−2 . Alayachi et al. [9] obtained the solutions of the following systems xn+1 = xn−3yn−4 yn(1 + xn−1yn−2xn−3yn−4) , yn+1 = yn−3xn−4 xn(±1 ±yn−1xn−2yn−3xn−4) . More information about systems of difference equations and their stability can be obtained in [10–16]. The motivation of this paper arises from the investigation of system (1.1). Hence, the main purpose of this paper is to investigate the solutions of the following system using Padovan numbers: Φn+1 = Φn−3 + 1 Ψn−1Φn−3 , Ψn+1 = Ψn−3 + 1 Φn−1Ψn−3 , (1.2) Φn+1 = Φn−3 − 1 Ψn−1Φn−3 , Ψn+1 = Ψn−3 − 1 Φn−1Ψn−3 , (1.3) where the initial conditions are real numbers. We also present the future pattern of the considered systems under some selected values for the initial conditions. The Padovan sequence {ρn}n∈N is defined by ρn+1 = ρn−1 + ρn−2, where ρ−2 = 0, ρ−1 = 0, ρ0 = 1. Int. J. Anal. Appl. (2022), 20:38 3 This paper is organized as follows. Section 4 presents the solutions of system (1.2) while Section 3 gives the solutions of system (1.3). Section (4) analyzes the stability of the equilibrium points. The numerical solutions of the considered systems are shown in Section 5. Finally, we conclude this article in Section 6. 2. Solutions of System (1.2) This section is devoted to give explicit forms for the solutions of system (1.2) using Padovan numbers. Theorem 2.1. Let {Φn, Ψn}∞n=−3 be the solutions of system (1.2) and assume that Φ−3 = α, Φ−2 = β, Φ−1 = γ, Φ0 = δ, Ψ−3 = �, Ψ−2 = ζ, Ψ−1 = η, and Ψ0 = κ. Then, for n = 0, 1, ..., we have Φ4n−3 = α (ρ2n−1η + ρ2n) + ρ2n−2 α (ρ2n−2η + ρ2n−1) + ρ2n−3 , Φ4n−2 = β(ρ2n−1κ + ρ2n) + ρ2n−2 β(ρ2n−2κ + ρ2n−1) + ρ2n−3 , Φ4n−1 = �(ρ2nγ + ρ2n+1) + ρ2n−1 �(ρ2n−1γ + ρ2n) + ρ2n−2 , Φ4n = ζ (ρ2nδ + ρ2n+1) + ρ2n−1 ζ (ρ2n−1δ + ρ2n) + ρ2n−2 . And Ψ4n−3 = � (ρ2n−1γ + ρ2n) + ρ2n−2 � (ρ2n−2γ + ρ2n−1) + ρ2n−3 , Ψ4n−2 = ζ (ρ2n−1δ + ρ2n) + ρ2n−2 ζ (ρ2n−2δ + ρ2n−1) + ρ2n−3 , Ψ4n−1 = α(ρ2nη + ρ2n+1) + ρ2n−1 α(ρ2n−1η + ρ2n) + ρ2n−2 , Ψ4n = β(ρ2nκ + ρ2n+1) + ρ2n−1 β(ρ2n−1κ + ρ2n) + ρ2n−2 , where ρn+1 = ρn−1 + ρn−2 ∀n ∈N is the Padovan sequence. Proof. The forms of the solutions are true for n = 0. We assume that n > 0 and that our assumption holds for n− 1. That is, Φ4n−7 = α (ρ2n−3η + ρ2n−2) + ρ2n−4 α (ρ2n−4η + ρ2n−3) + ρ2n−5 , Φ4n−6 = β(ρ2n−3κ + ρ2n−2) + ρ2n−4 β(ρ2n−4κ + ρ2n−3) + ρ2n−5 , Φ4n−5 = �(ρ2n−2γ + ρ2n−1) + ρ2n−3 �(ρ2n−3γ + ρ2n−2) + ρ2n−4 , 4 Int. J. Anal. Appl. (2022), 20:38 Φ4n−4 = ζ(ρ2n−2δ + ρ2n−1) + ρ2n−3 ζ(ρ2n−3δ + ρ2n−2) + ρ2n−4 , and Ψ4n−7 = � (ρ2n−3γ + ρ2n−2) + ρ2n−4 � (ρ2n−4γ + ρ2n−3) + ρ2n−5 , Ψ4n−6 = ζ (ρ2n−3δ + ρ2n−2) + ρ2n−4 ζ (ρ2n−4δ + ρ2n−3) + ρ2n−5 , Ψ4n−5 = α(ρ2n−2η + ρ2n−1) + ρ2n−3 α(ρ2n−3η + ρ2n−2) + ρ2n−4 , Ψ4n−4 = β(ρ2n−2κ + ρ2n−1) + ρ2n−3 β(ρ2n−3κ + ρ2n−2) + ρ2n−4 . Next, it can be easily observed from system (1.2) that Φ4n−3 = Φ4n−7 + 1 Ψ4n−5Φ4n−7 = α (ρ2n−3η+ρ2n−2)+ρ2n−4 α (ρ2n−4η+ρ2n−3)+ρ2n−5 + 1( α(ρ2n−2η+ρ2n−1)+ρ2n−3 α(ρ2n−3η+ρ2n−2)+ρ2n−4 )( α (ρ2n−3η+ρ2n−2)+ρ2n−4 α (ρ2n−4η+ρ2n−3)+ρ2n−5 ) = α (ρ2n−3η+ρ2n−2)+ρ2n−4+α (ρ2n−4η+ρ2n−3)+ρ2n−5 α (ρ2n−4η+ρ2n−3)+ρ2n−5 α(ρ2n−2η+ρ2n−1)+ρ2n−3 α (ρ2n−4η+ρ2n−3)+ρ2n−5 = α (ρ2n−3η + ρ2n−2) + ρ2n−4 + α (ρ2n−4η + ρ2n−3) + ρ2n−5 α(ρ2n−2η + ρ2n−1) + ρ2n−3 = α ρ2n−3η + αρ2n−2 + ρ2n−4 + α ρ2n−4η + αρ2n−3 + ρ2n−5 α(ρ2n−2η + ρ2n−1) + ρ2n−3 = α η (ρ2n−3 + ρ2n−4) + α (ρ2n−2 + ρ2n−3) + (ρ2n−4 + ρ2n−5) α(ρ2n−2η + ρ2n−1) + ρ2n−3 . Since ρn+1 = ρn−1 + ρn−2, we have Φ4n−3 = α (ρ2n−1η + ρ2n) + ρ2n−2 α(ρ2n−2η + ρ2n−1) + ρ2n−3 . Ψ4n−3 = Ψ4n−7 + 1 Φ4n−5Ψ4n−7 = � (ρ2n−3γ+ρ2n−2)+ρ2n−4 � (ρ2n−4γ+ρ2n−3)+ρ2n−5 + 1( �(ρ2n−2γ+ρ2n−1)+ρ2n−3 �(ρ2n−3γ+ρ2n−2)+ρ2n−4 )( � (ρ2n−3γ+ρ2n−2)+ρ2n−4 � (ρ2n−4γ+ρ2n−3)+ρ2n−5 ) = � (ρ2n−3γ+ρ2n−2)+ρ2n−4+� (ρ2n−4γ+ρ2n−3)+ρ2n−5 � (ρ2n−4γ+ρ2n−3)+ρ2n−5 �(ρ2n−2γ+ρ2n−1)+ρ2n−3 � (ρ2n−4γ+ρ2n−3)+ρ2n−5 = � (ρ2n−3γ + ρ2n−2) + ρ2n−4 + � (ρ2n−4γ + ρ2n−3) + ρ2n−5 �(ρ2n−2γ + ρ2n−1) + ρ2n−3 = �ρ2n−3γ + �ρ2n−2 + ρ2n−4 + �ρ2n−4γ + �ρ2n−3 + ρ2n−5 �(ρ2n−2γ + ρ2n−1) + ρ2n−3 Int. J. Anal. Appl. (2022), 20:38 5 = �γ (ρ2n−3 + ρ2n−4) + � (ρ2n−2 + ρ2n−3) + ρ2n−4 + ρ2n−5 �(ρ2n−2γ + ρ2n−1) + ρ2n−3 = �γρ2n−1 + �ρ2n + ρ2n−2 �(ρ2n−2γ + ρ2n−1) + ρ2n−3 = �(ρ2n−1γ + ρ2n) + ρ2n−2 �(ρ2n−2γ + ρ2n−1) + ρ2n−3 . Also, system (1.2) gives Φ4n−2 = Φ4n−6 + 1 Ψ4n−4Φ4n−6 = β(ρ2n−3κ+ρ2n−2)+ρ2n−4 β(ρ2n−4κ+ρ2n−3)+ρ2n−5 + 1( β(ρ2n−2κ+ρ2n−1)+ρ2n−3 β(ρ2n−3κ+ρ2n−2)+ρ2n−4 )( β(ρ2n−3κ+ρ2n−2)+ρ2n−4 β(ρ2n−4κ+ρ2n−3)+ρ2n−5 ) = β(ρ2n−3κ+ρ2n−2)+ρ2n−4+β(ρ2n−4κ+ρ2n−3)+ρ2n−5 β(ρ2n−4κ+ρ2n−3)+ρ2n−5 β(ρ2n−2κ+ρ2n−1)+ρ2n−3 β(ρ2n−4κ+ρ2n−3)+ρ2n−5 = β(ρ2n−3κ + ρ2n−2) + ρ2n−4 + β(ρ2n−4κ + ρ2n−3) + ρ2n−5 β(ρ2n−2κ + ρ2n−1) + ρ2n−3 = βρ2n−3κ + βρ2n−2 + ρ2n−4 + βρ2n−4κ + βρ2n−3 + ρ2n−5 β(ρ2n−2κ + ρ2n−1) + ρ2n−3 = βκ (ρ2n−3 + ρ2n−4) + β (ρ2n−2 + ρ2n−3) + ρ2n−4 + ρ2n−5 β(ρ2n−2κ + ρ2n−1) + ρ2n−3 = βκρ2n−1 + βρ2n + ρ2n−2 β(ρ2n−2κ + ρ2n−1) + ρ2n−3 = β (ρ2n−1κ + ρ2n) + ρ2n−2 β(ρ2n−2κ + ρ2n−1) + ρ2n−3 . Ψ4n−2 = Ψ4n−6 + 1 Φ4n−4Ψ4n−6 = ζ (ρ2n−3δ+ρ2n−2)+ρ2n−4 ζ (ρ2n−4δ+ρ2n−3)+ρ2n−5 + 1( ζ(ρ2n−2δ+ρ2n−1)+ρ2n−3 ζ(ρ2n−3δ+ρ2n−2)+ρ2n−4 )( ζ (ρ2n−3δ+ρ2n−2)+ρ2n−4 ζ (ρ2n−4δ+ρ2n−3)+ρ2n−5 ) = ζ (ρ2n−3δ+ρ2n−2)+ρ2n−4+ζ (ρ2n−4δ+ρ2n−3)+ρ2n−5 ζ (ρ2n−4δ+ρ2n−3)+ρ2n−5 ζ(ρ2n−2δ+ρ2n−1)+ρ2n−3 ζ (ρ2n−4δ+ρ2n−3)+ρ2n−5 = ζ (ρ2n−3δ + ρ2n−2) + ρ2n−4 + ζ (ρ2n−4δ + ρ2n−3) + ρ2n−5 ζ(ρ2n−2δ + ρ2n−1) + ρ2n−3 = ζ ρ2n−3δ + ζ ρ2n−2 + ρ2n−4 + ζ ρ2n−4δ + ζ ρ2n−3 + ρ2n−5 ζ(ρ2n−2δ + ρ2n−1) + ρ2n−3 = ζδ(ρ2n−3 + ρ2n−4) + ζ (ρ2n−2 + ρ2n−3) + ρ2n−4 + ρ2n−5 ζ(ρ2n−2δ + ρ2n−1) + ρ2n−3 = ζδρ2n−1 + ζρ2n + ρ2n−2 ζ(ρ2n−2δ + ρ2n−1) + ρ2n−3 = ζ (ρ2n−1δ + ρ2n) + ρ2n−2 ζ(ρ2n−2δ + ρ2n−1) + ρ2n−3 . Furthermore, from system (1.2), we have 6 Int. J. Anal. Appl. (2022), 20:38 Φ4n−1 = Φ4n−5 + 1 Ψ4n−3Φ4n−5 = �(ρ2n−2γ+ρ2n−1)+ρ2n−3 �(ρ2n−3γ+ρ2n−2)+ρ2n−4 + 1( � (ρ2n−1γ+ρ2n)+ρ2n−2 � (ρ2n−2γ+ρ2n−1)+ρ2n−3 )( �(ρ2n−2γ+ρ2n−1)+ρ2n−3 �(ρ2n−3γ+ρ2n−2)+ρ2n−4 ) = �(ρ2n−2γ+ρ2n−1)+ρ2n−3+�(ρ2n−3γ+ρ2n−2)+ρ2n−4 �(ρ2n−3γ+ρ2n−2)+ρ2n−4 � (ρ2n−1γ+ρ2n)+ρ2n−2 �(ρ2n−3γ+ρ2n−2)+ρ2n−4 = �(ρ2n−2γ + ρ2n−1) + ρ2n−3 + �(ρ2n−3γ + ρ2n−2) + ρ2n−4 � (ρ2n−1γ + ρ2n) + ρ2n−2 = �(ρ2nγ + ρ2n+1) + ρ2n−1 � (ρ2n−1γ + ρ2n) + ρ2n−2 . Ψ4n−1 = Ψ4n−5 + 1 Φ4n−3Ψ4n−5 = α(ρ2n−2η+ρ2n−1)+ρ2n−3 α(ρ2n−3η+ρ2n−2)+ρ2n−4 + 1( α (ρ2n−1η+ρ2n)+ρ2n−2 α (ρ2n−2η+ρ2n−1)+ρ2n−3 )( α(ρ2n−2η+ρ2n−1)+ρ2n−3 α(ρ2n−3η+ρ2n−2)+ρ2n−4 ) = α(ρ2n−2η+ρ2n−1)+ρ2n−3+α(ρ2n−3η+ρ2n−2)+ρ2n−4 α(ρ2n−3η+ρ2n−2)+ρ2n−4 α (ρ2n−1η+ρ2n)+ρ2n−2 α(ρ2n−3η+ρ2n−2)+ρ2n−4 = α(ρ2nη + ρ2n+1) + ρ2n−1 α (ρ2n−1η + ρ2n) + ρ2n−2 . Moreover, from system (1.2), we have Φ4n = 4n−4 + 1 4n−2Φ4n−4 = ζ(ρ2n−2δ+ρ2n−1)+ρ2n−3 ζ(ρ2n−3δ+ρ2n−2)+ρ2n−4 + 1( ζ (ρ2n−1δ+ρ2n)+ρ2n−2 ζ (ρ2n−2δ+ρ2n−1)+ρ2n−3 )( ζ(ρ2n−2δ+ρ2n−1)+ρ2n−3 ζ(ρ2n−3δ+ρ2n−2)+ρ2n−4 ) = ζ(ρ2n−2δ+ρ2n−1)+ρ2n−3+ζ(ρ2n−3δ+ρ2n−2)+ρ2n−4 �(ρ2n−3γ+ρ2n−2)+ρ2n−4 ζ (ρ2n−1δ+ρ2n)+ρ2n−2 ζ(ρ2n−3δ+ρ2n−2)+ρ2n−4 = ζ(ρ2nδ + ρ2n+1) + ρ2n−1 ζ (ρ2n−1δ + ρ2n) + ρ2n−2 . Ψ4n = Ψ4n−4 + 1 Φ4n−2Ψ4n−4 = β(ρ2n−2κ+ρ2n−1)+ρ2n−3 β(ρ2n−3κ+ρ2n−2)+ρ2n−4 + 1( β(ρ2n−1κ+ρ2n)+ρ2n−2 β(ρ2n−2κ+ρ2n−1)+ρ2n−3 )( β(ρ2n−2κ+ρ2n−1)+ρ2n−3 β(ρ2n−3κ+ρ2n−2)+ρ2n−4 ) Int. J. Anal. Appl. (2022), 20:38 7 = β(ρ2n−2κ+ρ2n−1)+ρ2n−3+β(ρ2n−3κ+ρ2n−2)+ρ2n−4 α(ρ2n−3η+ρ2n−2)+ρ2n−4 β(ρ2n−1κ+ρ2n)+ρ2n−2 β(ρ2n−3κ+ρ2n−2)+ρ2n−4 = β(ρ2nκ + ρ2n+1) + ρ2n−1 β(ρ2n−1κ + ρ2n) + ρ2n−2 . The proof has been completed. 3. Solutions of System (1.3) In this section, we show analytic solutions for system (1.3) using Padovan numbers. Theorem 3.1. Let {Φn, Ψn}∞n=−3 be the solutions of system (1.3) and assume that Φ−3 = α, Φ−2 = β, Φ−1 = γ, Φ0 = δ, Ψ−3 = �, Ψ−2 = ζ, Ψ−1 = η, and Ψ0 = κ. Then, for n = 0, 1, ..., we have Φ4n−3 = α (ρ2n −ρ2n−1η) −ρ2n−2 α (ρ2n−2η −ρ2n−1) + ρ2n−3 , Φ4n−2 = β(ρ2n −ρ2n−1κ) −ρ2n−2 β(ρ2n−2κ−ρ2n−1) + ρ2n−3 , Φ4n−1 = �(ρ2n+1 −ρ2nγ) −ρ2n−1 �(ρ2n−1γ −ρ2n) + ρ2n−2 , Φ4n = ζ (ρ2n+1 −ρ2nδ) −ρ2n−1 ζ (ρ2n−1δ −ρ2n) + ρ2n−2 . And Ψ4n−3 = � (ρ2n −ρ2n−1γ) −ρ2n−2 � (ρ2n−2γ −ρ2n−1) + ρ2n−3 , Ψ4n−2 = ζ (ρ2n −ρ2n−1δ) −ρ2n−2 ζ (ρ2n−2δ −ρ2n−1) + ρ2n−3 , Ψ4n−1 = α(ρ2n+1 −ρ2nη) −ρ2n−1 α(ρ2n−1η −ρ2n) + ρ2n−2 , Ψ4n = β(ρ2n+1 −ρ2nκ) −ρ2n−1 β(ρ2n−1κ−ρ2n) + ρ2n−2 , where ρn+1 = ρn−1 + ρn−2 ∀n ∈N is the Padovan sequence. Proof. The above solutions are true for n = 0. We suppose that n > 0 and that our assumption holds for n− 1. That is, Φ4n−7 = α (ρ2n−2 −ρ2n−3η) −ρ2n−4 α (ρ2n−4η −ρ2n−3) + ρ2n−5 , Φ4n−6 = β(ρ2n−2 −ρ2n−3κ) −ρ2n−4 β(ρ2n−4κ−ρ2n−3) + ρ2n−5 , Φ4n−5 = �(ρ2n−1 −ρ2n−2γ) −ρ2n−3 �(ρ2n−3γ −ρ2n−2) + ρ2n−4 , Φ4n−4 = ζ(ρ2n−1 −ρ2n−2δ) −ρ2n−3 ζ(ρ2n−3δ −ρ2n−2) + ρ2n−4 . 8 Int. J. Anal. Appl. (2022), 20:38 And Ψ4n−7 = � (ρ2n−2 −ρ2n−3γ) −ρ2n−4 � (ρ2n−4γ −ρ2n−3) + ρ2n−5 , Ψ4n−6 = ζ (ρ2n−2 −ρ2n−3δ) −ρ2n−4 ζ (ρ2n−4δ −ρ2n−3) + ρ2n−5 , Ψ4n−5 = α(ρ2n−1 −ρ2n−2η) −ρ2n−3 α(ρ2n−3η −ρ2n−2) + ρ2n−4 , Ψ4n−4 = β(ρ2n−1 −ρ2n−2κ) −ρ2n−3 β(ρ2n−3κ−ρ2n−2) + ρ2n−4 . System (1.3) leads to Φ4n−3 = Φ4n−7 − 1 Ψ4n−5Φ4n−7 = α (ρ2n−2−ρ2n−3η)−ρ2n−4 α (ρ2n−4η−ρ2n−3)+ρ2n−5 − 1( α(ρ2n−1−ρ2n−2η)−ρ2n−3 α(ρ2n−3η−ρ2n−2)+ρ2n−4 )( α(ρ2n−2−ρ2n−3η)−ρ2n−4 α (ρ2n−4η−ρ2n−3)+ρ2n−5 ) = α (ρ2n−2−ρ2n−3η)−ρ2n−4−(α (ρ2n−4η−ρ2n−3)+ρ2n−5) α (ρ2n−4η−ρ2n−3)+ρ2n−5 − α(ρ2n−1−ρ2n−2η)−ρ2n−3 α (ρ2n−4η−ρ2n−3)+ρ2n−5 = α (ρ2n−2 −ρ2n−3η) −ρ2n−4 − (α (ρ2n−4η −ρ2n−3) + ρ2n−5) −α(ρ2n−1 −ρ2n−2η) −ρ2n−3 = α (ρ2n−2 + ρ2n−3) −αη (ρ2n−3 + ρ2n−4) − (ρ2n−4 + ρ2n−5) α(ρ2n−2η −ρ2n−1) + ρ2n−3 = αρ2n −αηρ2n−1 −ρ2n−2 α(ρ2n−2η −ρ2n−1) + ρ2n−3 = α (ρ2n −ρ2n−1η) −ρ2n−2 α(ρ2n−2η −ρ2n−1) + ρ2n−3 . Ψ4n−3 = Ψ4n−7 − 1 Φ4n−5Ψ4n−7 = � (ρ2n−2−ρ2n−3γ)−ρ2n−4 � (ρ2n−4γ−ρ2n−3)+ρ2n−5 − 1( �(ρ2n−1−ρ2n−2γ)−ρ2n−3 �(ρ2n−3γ−ρ2n−2)+ρ2n−4 )( � (ρ2n−2−ρ2n−3γ)−ρ2n−4 � (ρ2n−4γ−ρ2n−3)+ρ2n−5 ) = � (ρ2n−2−ρ2n−3γ)−ρ2n−4−(� (ρ2n−4γ−ρ2n−3)+ρ2n−5) � (ρ2n−4γ−ρ2n−3)+ρ2n−5 − �(ρ2n−1−ρ2n−2γ)−ρ2n−3 � (ρ2n−4γ−ρ2n−3)+ρ2n−5 = � (ρ2n−2 −ρ2n−3γ) −ρ2n−4 − (� (ρ2n−4γ −ρ2n−3) + ρ2n−5) −�(ρ2n−1 −ρ2n−2γ) −ρ2n−3 = � (ρ2n−2 + ρ2n−3) − �γ (ρ2n−3 + ρ2n−4) − (ρ2n−4 + ρ2n−5) �(ρ2n−2γ −ρ2n−1) + ρ2n−3 = �ρ2n − �γρ2n−1 −ρ2n−2 �(ρ2n−2γ −ρ2n−1) + ρ2n−3 = � (ρ2n −ρ2n−1γ) −ρ2n−2 �(ρ2n−2γ −ρ2n−1) + ρ2n−3 . Furthermore, we can obtain from system (1.3) that Int. J. Anal. Appl. (2022), 20:38 9 Φ4n−2 = Φ4n−6 − 1 Ψ4n−4Φ4n−6 = β(ρ2n−2−ρ2n−3κ)−ρ2n−4 β(ρ2n−4κ−ρ2n−3)+ρ2n−5 − 1( β(ρ2n−1−ρ2n−2κ)−ρ2n−3 β(ρ2n−3κ−ρ2n−2)+ρ2n−4 )( β(ρ2n−2−ρ2n−3κ)−ρ2n−4 β(ρ2n−4κ−ρ2n−3)+ρ2n−5 ) = β(ρ2n−2−ρ2n−3κ)−ρ2n−4−(β(ρ2n−4κ−ρ2n−3)+ρ2n−5) β(ρ2n−4κ−ρ2n−3)+ρ2n−5 −β(ρ2n−1−ρ2n−2κ)−ρ2n−3 β(ρ2n−4κ−ρ2n−3)+ρ2n−5 = β(ρ2n−2 −ρ2n−3κ) −ρ2n−4 − (β(ρ2n−4κ−ρ2n−3) + ρ2n−5) −β(ρ2n−1 −ρ2n−2κ) −ρ2n−3 = β(ρ2n −ρ2n−1κ) −ρ2n−2 β(ρ2n−2κ−ρ2n−1) + ρ2n−3 . Ψ4n−2 = Ψ4n−6 − 1 Φ4n−4Ψ4n−6 = ζ (ρ2n−2−ρ2n−3δ)−ρ2n−4 ζ (ρ2n−4δ−ρ2n−3)+ρ2n−5 − 1( ζ(ρ2n−1−ρ2n−2δ)−ρ2n−3 ζ(ρ2n−3δ−ρ2n−2)+ρ2n−4 )( ζ (ρ2n−2−ρ2n−3δ)−ρ2n−4 ζ (ρ2n−4δ−ρ2n−3)+ρ2n−5 ) = ζ (ρ2n−2−ρ2n−3δ)−ρ2n−4−(ζ (ρ2n−4δ−ρ2n−3)+ρ2n−5) ζ (ρ2n−4δ−ρ2n−3)+ρ2n−5 − ζ(ρ2n−1−ρ2n−2δ)−ρ2n−3 ζ (ρ2n−4δ−ρ2n−3)+ρ2n−5 = ζ (ρ2n−2 −ρ2n−3δ) −ρ2n−4 − (ζ (ρ2n−4δ −ρ2n−3) + ρ2n−5) −ζ(ρ2n−1 −ρ2n−2δ) −ρ2n−3 = ζ (ρ2n −ρ2n−1δ) −ρ2n−2 ζ(ρ2n−2δ −ρ2n−1) + ρ2n−3 . In addition, one can obtain from system (1.3) that Φ4n−1 = Φ4n−5 − 1 Ψ4n−3Φ4n−5 = �(ρ2n−1−ρ2n−2γ)−ρ2n−3 �(ρ2n−3γ−ρ2n−2)+ρ2n−4 − 1( � (ρ2n−ρ2n−1γ)−ρ2n−2 � (ρ2n−2γ−ρ2n−1)+ρ2n−3 )( �(ρ2n−1−ρ2n−2γ)−ρ2n−3 �(ρ2n−3γ−ρ2n−2)+ρ2n−4 ) = �(ρ2n−1−ρ2n−2γ)−ρ2n−3−(�(ρ2n−3γ−ρ2n−2)+ρ2n−4) �(ρ2n−3γ−ρ2n−2)+ρ2n−4 − � (ρ2n−ρ2n−1γ)−ρ2n−2 �(ρ2n−3γ−ρ2n−2)+ρ2n−4 = �(ρ2n−1 −ρ2n−2γ) −ρ2n−3 − (�(ρ2n−3γ −ρ2n−2) + ρ2n−4) −� (ρ2n −ρ2n−1γ) −ρ2n−2 = �(ρ2n+1 −ρ2nγ) −ρ2n−1 � (ρ2n−1γ −ρ2n) + ρ2n−2 . 10 Int. J. Anal. Appl. (2022), 20:38 Ψ4n−1 = Ψ4n−5 − 1 Φ4n−3Ψ4n−5 = α(ρ2n−1−ρ2n−2η)−ρ2n−3 α(ρ2n−3η−ρ2n−2)+ρ2n−4 − 1( α (ρ2n−ρ2n−1η)−ρ2n−2 α (ρ2n−2η−ρ2n−1)+ρ2n−3 )( α(ρ2n−1−ρ2n−2η)−ρ2n−3 α(ρ2n−3η−ρ2n−2)+ρ2n−4 ) = α(ρ2n−1−ρ2n−2η)−ρ2n−3−(α(ρ2n−3η−ρ2n−2)+ρ2n−4) α(ρ2n−3η−ρ2n−2)+ρ2n−4 − α (ρ2n−ρ2n−1η)−ρ2n−2 α(ρ2n−3η−ρ2n−2)+ρ2n−4 = α(ρ2n−1 −ρ2n−2η) −ρ2n−3 − (α(ρ2n−3η −ρ2n−2) + ρ2n−4) −α (ρ2n −ρ2n−1η) −ρ2n−2 = α(ρ2n+1 −ρ2nη) −ρ2n−1 α (ρ2n−1η −ρ2n) + ρ2n−2 . Finally, system (1.3) leads to Φ4n = Φ4n−4 − 1 Ψ4n−2Φ4n−4 = ζ(ρ2n−1−ρ2n−2δ)−ρ2n−3 ζ(ρ2n−3δ−ρ2n−2)+ρ2n−4 − 1( ζ (ρ2n−ρ2n−1δ)−ρ2n−2 ζ (ρ2n−2δ−ρ2n−1)+ρ2n−3 )( ζ(ρ2n−1−ρ2n−2δ)−ρ2n−3 ζ(ρ2n−3δ−ρ2n−2)+ρ2n−4 ) = ζ(ρ2n−1−ρ2n−2δ)−ρ2n−3−(ζ(ρ2n−3δ−ρ2n−2)+ρ2n−4) ζ(ρ2n−3δ−ρ2n−2)+ρ2n−4 − ζ (ρ2n−ρ2n−1δ)−ρ2n−2 ζ(ρ2n−3δ−ρ2n−2)+ρ2n−4 = ζ(ρ2n−1 −ρ2n−2δ) −ρ2n−3 − (ζ(ρ2n−3δ −ρ2n−2) + ρ2n−4) −ζ (ρ2n −ρ2n−1δ) −ρ2n−2 = ζ(ρ2n+1 −ρ2nδ) −ρ2n−1 ζ (ρ2n−1δρ2n) + ρ2n−2 . Ψ4n = Ψ4n−4 − 1 Φ4n−2Ψ4n−4 = β(ρ2n−1−ρ2n−2κ)−ρ2n−3 β(ρ2n−3κ−ρ2n−2)+ρ2n−4 − 1( β(ρ2n−ρ2n−1κ)−ρ2n−2 β(ρ2n−2κ−ρ2n−1)+ρ2n−3 )( β(ρ2n−1−ρ2n−2κ)−ρ2n−3 β(ρ2n−3κ−ρ2n−2)+ρ2n−4 ) = β(ρ2n−1−ρ2n−2κ)−ρ2n−3−(β(ρ2n−3κ−ρ2n−2)+ρ2n−4) β(ρ2n−3κ−ρ2n−2)+ρ2n−4 − β(ρ2n−ρ2n−1κ)−ρ2n−2 β(ρ2n−3κ−ρ2n−2)+ρ2n−4 = β(ρ2n−1 −ρ2n−2κ) −ρ2n−3 − (β(ρ2n−3κ−ρ2n−2) + ρ2n−4) −β(ρ2n −ρ2n−1κ) −ρ2n−2 = β(ρ2n+1 −ρ2nκ) −ρ2n−1 β(ρ2n−1κ−ρ2n) + ρ2n−2 . Hence, the proof is completed. Int. J. Anal. Appl. (2022), 20:38 11 4. Stability of the equilibrium points This section discusses the stability of the equilibrium points of the considered systems. Theorem 4.1. System (1.2) has a unique real equilibrium point (q1,q1) which is a saddle point. Proof. The equilibrium point of system (1.2) is given by Ψ∗ = Ψ∗ + 1 Φ∗Ψ∗ , Φ∗ = Φ∗ + 1 Φ∗Ψ∗ . (4.1) Subtract the second equation of system (4.1) from the first equation to have Φ∗ 2 Ψ∗ − Φ∗ − 1 − Φ∗ Ψ∗ 2 + Ψ∗ + 1 = 0, or, (Φ∗ Ψ∗ − 1) (Φ∗ − Ψ∗) = 0. Hence, Φ∗ Ψ∗ − 1 = 0 , Φ∗ Ψ∗ = 1 , (4.2) Φ∗ − Ψ∗ = 0 , Ψ∗ = Φ∗ . (4.3) Equation. (4.2) dose not satisfy system (4.1). Therefore, the only real equilibrium point is obtained from substituting equation. (4.3) into any equation of system (4.1). This gives Φ∗ 3 − Φ∗ − 1 = 0. (4.4) Solving equation. (4.4) gives q1 = α2 + 12 6α , q2,3 = − α2 + 12 6α ± √ 3 2 ( α 6 − 2 3α ) i, where, α = 3 √ 108 + 12 √ 69. Hence, the only real equilibrium point is (q1,q1). Now, we find the Jacobian matrix. Let F (u,v) = (f (u,v),g(u,v)), where f (u,v) = u + 1 uv , g(u,v) = v + 1 uv . Then, JF =   −1 u2v −u+1 uv2 −v+1 u2v −1 uv2   . 12 Int. J. Anal. Appl. (2022), 20:38 Evaluating the Jacobian matrix about (q1,q1) gives, JF (q1,q1) =   −1 q31 −q1+1 q31 −q1+1 q31 −1 q31   . Thus, the characteristic equation of this matrix is given by ( λ + 1 q31 )( λ + 1 q31 ) − (q1 + 1) 2 q61 = 0. Or, ( λ + 1 q31 )2 − (q1 + 1) 2 q61 = 0. Note that (q1+1) 2 q61 = 1. Hence, ( λ + 1 q31 )2 − 1 = 0. Then, λ + 1 q31 = ±1. For λ + 1 q31 = 1, we have λ1 = |1− 1q31 | < 1. For λ + 1 q31 = −1, we have λ2 = |−1− 1q31 | = 1 + 1 q31 > 1. Since λ1 < 1 and λ2 > 1, the point (q1,q1) is a saddle point. Theorem 4.2. System (1.3) has a unique real equilibrium point (−q1,−q1) which is a saddle point. Proof. The proof is similar to the proof of Theorem 4.1. 5. Behavior of the solutions This section presents the future pattern of the considered systems under specific initial conditions. We selected some random values for the initial conditions to illustrate the long behavior solutions. Example 1. Figure 1 (left) presents the dynamical behavior of the solutions of system (1.2) under the selected values Φ−3 = 7, Φ−2 = 3, Φ−1 = 4, Φ0 = 6, Ψ−3 = 5, Ψ−2 = 3, Ψ−1 = 5, and Ψ0 = 6. Example 2. The exact solutions of system (1.2) are also plotted in Figure 1 (right) under the initial conditions Φ−3 = 1, Φ−2 = 3, Φ−1 = 1, Φ0 = 2, Ψ−3 = 1, Ψ−2 = 4, Ψ−1 = 2, and Ψ0 = 4. Int. J. Anal. Appl. (2022), 20:38 13 0 10 20 30 40 50 n 0 1 2 3 4 5 6 7 Plot of The First System x(n) y(n) 0 10 20 30 40 50 n 0 0.5 1 1.5 2 2.5 3 3.5 4 Plot of The First System x(n) y(n) Figure 1. The dynamical behavior of the solutions of system (1.2). Example 3. In Figure 2 (left), we plot the dynamical behavior of the solutions of system (1.3) under the selected values Φ−3 = 7, Φ−2 = 2.5, Φ−1 = 6, Φ0 = 4.4, Ψ−3 = 3.8, Ψ−2 = 2, Ψ−1 = 5, and Ψ0 = 3.3. Example 4. This example illustrates the behavior of the exact solutions of system (1.3) when we use the following initial values: Φ−3 = 7, Φ−2 = 3, Φ−1 = 4, Φ0 = 6, Ψ−3 = 5, Ψ−2 = 3, Ψ−1 = 5, and Ψ0 = 6. See Figure 2 (right). 0 5 10 15 20 25 30 35 40 45 50 n -10 -8 -6 -4 -2 0 2 4 6 8 Plot of The Second System x(n) y(n) 0 5 10 15 20 25 30 35 40 45 50 n -4 -2 0 2 4 6 8 Plot of The Second System x(n) y(n) Figure 2. The dynamical behavior of the solutions of system (1.3). 6. Conclusion In brief, we have discussed the solutions of the considered systems using Padovan numbers. Succes- sive iterations have been successfully used in extracting the exact solutions. Theorem 2.1 presents the solutions of system 1.2 while Theorem 3.1 gives the solutions of system 1.3. The obtained solutions 14 Int. J. Anal. Appl. (2022), 20:38 are presented in the form of rational relations. Jacobian matrix has been successfully used to examine the stability of the real equilibrium points. The equilibrium points are saddle. In Section 5, we depict the behavior of the solutions for some random initial conditions. The constructed exact solutions are stable. The used approaches can be utilized to deal with high order systems of difference equations. Conflicts of Interest: The author declares that there are no conflicts of interest regarding the publi- cation of this paper. References [1] J.D. Murray, Mathematical Biology: I. An Introduction, 3rd Ed., Springer-Verlag, New York, 2001. https://doi.org/10.1007/b98868. [2] D.T. Tollu, Y. Yazlik, N. Taskara, On the Solutions of Two Special Types of Riccati Difference Equation via Fibonacci Numbers, Adv. Differ. Equ. 2013 (2013), 174. https://doi.org/10. 1186/1687-1847-2013-174. [3] Y. Yazlik, D.T. Tollu, N. Taskara, On the Solutions of Difference Equation Systems with Padovan Numbers, Appl. Math. 04 (2013), 15–20. https://doi.org/10.4236/am.2013.412a002. [4] M.B. Almatrafi, Solutions Structures for Some Systems of Fractional Difference Equations, Open J. Math. Anal. 3 (2019), 52–61. https://doi.org/10.30538/psrp-oma2019.0032. [5] M.B. Almatrafi, M.M. Alzubaidi, Analysis of the Qualitative Behaviour of an Eighth-Order Frac- tional Difference Equation, Open J. Discret. Appl. Math. 2 (2019), 41–47. https://doi.org/ 10.30538/psrp-odam2019.0010. [6] C. Çinar, On the Positive Solutions of the Difference Equation System xn+1 = 1 yn , yn+1 = yn yn−1xn−1 , Appl. Math. Comput. 158 (2004), 303–305. https://doi.org/10.1016/j.amc.2003. 08.073. [7] E.M. Elsayed, Solutions of Rational Difference Systems of Order Two, Math. Computer Model. 55 (2012), 378–384. https://doi.org/10.1016/j.mcm.2011.08.012. [8] H.S. Alayachi, A.Q. Khan, M.S.M. Noorani, On the Solutions of Three-Dimensional Rational Difference Equation Systems, J. Math. 2021 (2021), 2480294. https://doi.org/10.1155/ 2021/2480294. [9] H.S. Alayachi, A.Q. Khan, M.S.M. Noorani, A. Khaliq, Displaying the Structure of the Solutions for Some Fifth-Order Systems of Recursive Equations, Math. Probl. Eng. 2021 (2021), 6682009. https://doi.org/10.1155/2021/6682009. [10] H.S. Alayachi, M.S.M. Noorani, A.Q. Khan, M.B. Almatrafi, Analytic Solutions and Stability of Sixth Order Difference Equations, Math. Probl. Eng. 2020 (2020), 1230979. https://doi.org/ 10.1155/2020/1230979. [11] S. Elaydi, An Introduction to Difference Equations, Springer-Verlag, New York, 2005. https: //doi.org/10.1007/0-387-27602-5. https://doi.org/10.1007/b98868 https://doi.org/10.1186/1687-1847-2013-174 https://doi.org/10.1186/1687-1847-2013-174 https://doi.org/10.4236/am.2013.412a002 https://doi.org/10.30538/psrp-oma2019.0032 https://doi.org/10.30538/psrp-odam2019.0010 https://doi.org/10.30538/psrp-odam2019.0010 https://doi.org/10.1016/j.amc.2003.08.073 https://doi.org/10.1016/j.amc.2003.08.073 https://doi.org/10.1016/j.mcm.2011.08.012 https://doi.org/10.1155/2021/2480294 https://doi.org/10.1155/2021/2480294 https://doi.org/10.1155/2021/6682009 https://doi.org/10.1155/2020/1230979 https://doi.org/10.1155/2020/1230979 https://doi.org/10.1007/0-387-27602-5 https://doi.org/10.1007/0-387-27602-5 Int. J. Anal. Appl. (2022), 20:38 15 [12] M.B. Almatrafi, E.M. Elsayed, F. Alzahrani, Qualitative Behavior of Two Rational Difference Equations, Fund. J. Math. Appl. 1 (2018), 194–204. https://doi.org/10.33401/fujma. 454999. [13] M.B. Almatrafi, E.M. Elsayed, Solutions And Formulae For Some Systems Of Difference Equa- tions, MathLAB J. 1 (2018), 356-369. [14] M.B. Almatrafi, E.M. Elsayed, F. Alzahrani, Qualitative Behavior of a Quadratic Second-Order Rational Difference Equation, Int. J. Adv. Math. 2019 (2019), 1-14. [15] M.B. Almatrafi, Exact Solutions and Stability of Sixth Order Difference Equations, Electron. J. Math. Anal. Appl. 10 (2022), 209-225. [16] M.B. Almatrafi, Abundant Traveling Wave and Numerical Solutions for Novikov-Veselov System With Their Stability and Accuracy, Appl. Anal. (2022). https://doi.org/10.1080/00036811. 2022.2027381. https://doi.org/10.33401/fujma.454999 https://doi.org/10.33401/fujma.454999 https://doi.org/10.1080/00036811.2022.2027381 https://doi.org/10.1080/00036811.2022.2027381 1. Introduction 2. Solutions of System (1.2) 3. Solutions of System (1.3) 4. Stability of the equilibrium points 5. Behavior of the solutions 6. Conclusion References