EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6506 ISSN 1307-5543 – ejpam.com Published by New York Business Global On the Dynamics of Particular Discrete Quadratic Systems with Potential Absence of Equilibria: Standard Euler versus Nonstandard Mickens Discretizations Moch. Fandi Ansori1, Hafidh Khoerul Fata1,∗, Nurcahya Yulian Ashar1, Ratna Herdiana1, Khairul Saleh2 1 Department of Mathematics, Faculty of Science and Mathematics, Universitas Diponegoro, Semarang, 50275, Indonesia 2 Department of Mathematics, College of Computing and Mathematics, King Fahd University of Petroleum and Minerals, Dhahran, 31261, Saudi Arabia Abstract. This study examined the dynamics of a three-dimensional discrete-time quadratic sys- tem that may lack equilibrium points, depending on the parameter values. Two discretization methods were considered: the standard forward Euler method and a nonstandard finite-difference scheme following Mickens’ approach. We investigated how the control parameter and step size af- fect stability and transition to chaos. Analytical derivations, together with numerical experiments including phase portraits, bifurcation diagrams, and Lyapunov exponent calculations, reveal di- verse dynamic behaviors. Both discretizations capture chaotic motion, but the nonstandard finite- difference scheme preserves structural stability more effectively across parameter ranges. These results enrich the study of discrete chaotic systems and highlight the advantages of nonstandard discretization. 2020 Mathematics Subject Classifications: 37M15, 39A33 Key Words and Phrases: Discrete-time dynamical systems, quadratic maps, nonstandard finite difference scheme, two-parameter bifurcation analysis, Lyapunov exponent analysis 1. Introduction Dynamic systems, whether continuous or discrete, are used in many fields such as physics and engineering [1–4], chemistry and environmental science [5], biology and epi- demiology [6–9], and economics and finance [10–12]. Analyzing these systems is challeng- ing because chaotic behavior may arise even in low-dimensional models with only three differential equations [13]. A notable example is a three-dimensional dynamical system that may lack equilibria [14]. Such systems, characterized by the absence of fixed points, can be studied using ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6506 Email addresses: hafid@lecturer.undip.ac.id (H. K. Fata) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 2 of 27 both standard and non-standard difference schemes. These approaches provide insights into their complex behavior and help bridge mathematical analyses with potential real- world applications [15, 16]. Among the many chaotic systems, an especially intriguing class is those without equi- libria. Jafari [14] conducted a comprehensive investigation of various instances of three- dimensional continuous-time systems without equilibria, each exhibiting quadratic nonlin- earities in a relatively simple form. A total of 17 such systems were identified and labeled NE1 through NE17, where ”NE” denotes ”no equilibria.” In this study, we consider model NE11 in [14] for further discrete-time dynamic analysis, which is given by the following dynamical system  ẋ = y ẏ = −x+ z ż = z − 2xy − 1.8xz − a, (1) where (x0, y0, z0) = (0, 1.6, 3). Jafari [14] showed that model (1) with a = 1 exhibits chaotic behavior. Beyond its theoretical interest, the study of such systems is motivated by practical applications. Three-dimensional systems without equilibrium points, such as NE11, also have broad potential in modern technology applications. For example, Gong et al. [17] developed a chaos system with hidden attractors, which was applied to random-number generator-based image encryption, thereby proving its direct relevance to digital security. In the field of applied physics, Djorwe et al. [18] explored hidden attractors in optome- chanical systems, whereas Kuznetsov et al. [19] presented a theoretical and experimental approach for hidden attractors in Chua’s circuits. Furthermore, Yang and Lai [20] designed a discrete memristor-based hyperchaotic map and demonstrated its potential in electronic and neuromorphic designs. From a methodological perspective, Maaita and Prousalis [21] compared four chaos indicators in systems with hidden attractors, which are useful for the numerical validation of discrete models. Most recently, a study on multi-scroll hidden at- tractors [22] showed great potential in cryptography and information security. Thus, it is interesting to explore the NE11 model in a discrete context, not only from a mathematical point of view but also from a real-world application perspective. A key feature that links both theoretical analysis and technological applications is the notion of hidden attractors that are not linked to equilibrium [23]. These attractors highlight the richness of dynamics that can emerge without fixed points. Pham [2] also examined hidden attractors in discrete-time systems and showed that complex behavior may appear even when stable equilibria exist [24]. This phenomenon has been documented across various dimensions, ranging from three-dimensional systems [24] to more complex higher-dimensional hyperchaotic systems [25]. To capture such phenomena in discrete-time, appropriate discretization schemes are needed. In this study, two distinct methodologies are employed to transform the continuous- time model (1) into a discrete-time model. The first methodology is the standard difference scheme, commonly known as the Euler scheme. The second methodology is a nonstandard difference scheme, which preserves the same equilibria as the Euler method but modifies M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 3 of 27 the discretization to improve stability. Nonstandard difference schemes are pivotal in the numerical analysis of dynamic sys- tems. Recent research has focused on the development of efficient non-standard finite- difference schemes tailored to various classes of dynamic systems. For instance, Nonlao- pon et al. [26] proposed an NSFD approach for predator–prey models, whereas Din et al. [27] applied this method to both classical and fractional-order SEIVR epidemic models. Kaur and Singh [28] introduced a hybrid scheme combining NSFD with a compact finite difference method to solve the one-dimensional Burgers’ equation. Khan et al. [29] inves- tigated COVID-19 transmission dynamics using a nonstandard discretization framework. In a related study, Gümüş and Türk [30] analyzed the behavior of a hepatitis B epidemic model and developed an NSFD scheme to capture its dynamics. These schemes facilitate high-resolution numerical solutions by meticulously designing the discretization process, thereby enhancing the capture of chaotic behavior and obtaining solutions that preserve the essential characteristics of the original equations [31]. The nonstandard finite difference scheme introduced by Mickens [32] is particularly useful when the right-hand side of a differential equation includes terms involving the same variable as the left-hand side, but with a negative sign. In the third equation of the system (1), the term −1.8xz includes the variable z, which appears as ż on the left-hand side, thereby satisfying Mickens’ criterion. This allows the system (1) to be appropri- ately transformed using the nonstandard scheme. While the original system is defined in continuous time, this study focuses on the discrete map resulting from nonstandard discretization, not merely as a numerical approximation, but also as an autonomous dy- namical system with its own qualitative behavior. Discrete-time dynamics often differ significantly from continuous-time flow. Here, we treat the discrete map as a system of intrinsic interest, focusing on how discretization affects chaos, sensitivity to step size and parameters, and the structure of attractors. Although nonstandard schemes are widely used to preserve the key features of chaotic systems, their application to three-dimensional quadratic maps remains limited. This study addresses this gap by analyzing the discrete dynamics derived from system (1). The remainder of this paper is organized as follows. Section 2 introduces the discrete- time formulation of the original continuous system using a standard difference scheme followed by numerical simulations to illustrate its dynamic behavior. Section 3 presents the discretization using Mickens’ nonstandard difference scheme, and investigates its qual- itative features through simulations. A comparative analysis highlighted the effects of discretization choice on the dynamics of the system. Finally, Section 4 summarizes the key findings and outlines the potential directions for future research. 2. Discrete-Time System via Standard Euler Scheme We discretize system (1) using the forward Euler scheme, where ẋ = xt+1−xt h for h > 0, yielding M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 4 of 27  xt+1 = xt + hyt yt+1 = yt + h(−xt + zt) zt+1 = zt + h(zt − 2xtyt − 1.8xtzt − a). (2) 2.1. Equilibrium Analysis and Stability Conditions The system (2) has two equilibria, E∗ = (x∗, y∗, z∗), namely E∗ 1 = ( 1 + √ 1− 7.2a 3.6 , 0, 1 + √ 1− 7.2a 3.6 ) , E∗ 2 = ( 1− √ 1− 7.2a 3.6 , 0, 1− √ 1− 7.2a 3.6 ) . If a > 1 7.2 , then the system (2) has no equilibrium. Therefore, in [14], model (1), with a = 1 ≥ 1 7.2 has no equilibrium. However, in this study, we examined a general case, that is, an arbitrary a ∈ R. Note that the equilibrium of variable y is always zero and z∗ = x∗. So, we evaluate the Jacobian matrix of system (2) at any point E = (x, 0, x) ∈ R3 as follows J(E) =  1 h 0 −h 1 h −1.8hx −2hx 1 + h− 1.8hx  . (3) The characteristic equation of J(E) is given by p(λ) = λ3 + b1λ 2 + b2λ+ b3, (4) where b1 = −(3 + h) + 1.8hx, b2 = (3 + 2h+ h2) + (−3.6h+ 2h2)x, b3 = −(1 + h+ h2 + h3) + (1.8h− 2h2 + 3.6h3)x. Following Andaluz [33], the Jury stability conditions for a three-dimensional discrete system require that an equilibrium point E = (x, 0, x) ∈ R3 is locally asymptotically stable if  1 + b1 + b2 + b3 > 0, 1− b1 + b2 − b3 > 0, 1− b2 + b1b3 − b23 > 0, b2 > 3. (5) First, we study equilibrium E∗ 2 . From 1 + b1 + b2 + b3 = h3(3.6x − 1) > 0, we have x > 1 3.6 . However, the value x ∗ = 1− √ 1−7.2a 3.6 in E∗ 2 was always less than 1 3.6 . Thus, E ∗ 2 fails to satisfy Jury’s conditions, and is always unstable. For the case of equilibrium E∗ 1 , we have the following theorem. M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 5 of 27 Theorem 1. Equilibrium E∗ 1 = (x∗, 0, x∗), where x∗ = 1 + √ 1− 7.2a 3.6 , is locally asymptotically stable if the following conditions are satisfied: 1 7.2 [ 1− ( − 1 + 3.6min{ 2+h 3.6+2h , (h+2)(h2+4) 7.2h−4h2+3.6h3 } )2] < a < 1 7.2 , c1 ( 1 + √ 1− 7.2a 3.6 )2 + c2 ( 1 + √ 1− 7.2a 3.6 ) + c3 > 0, (6) where  c1 = −12.96h6 + 14.4h5 − 10.48h4 + 3.6h3, c2 = 7.2h6 + 3.2h5 + 1.4h4 − 3.8h3, c3 = −h4 − 2h5. (7) Proof. The condition 1 + b1 + b2 + b3 = h3(3.6x∗ − 1) > 0 or x∗ > 1 3.6 is fulfilled. This condition is equivalent to a < 1 7.2 . (8) From the condition b2 > 3, we have x∗ < 2+h 3.6+2h and from 1− b1 + b2 − b3 = − ( 7.2h− 4h2 + 3.6h3 ) x∗ + (h+ 2)(h2 + 4) > 0, we have x∗ < (h+2)(h2+4) 7.2h−4h2+3.6h3 . These two conditions are equivalent to x∗ < min { 2 + h 3.6 + 2h , (h+ 2)(h2 + 4) 7.2h− 4h2 + 3.6h3 } or a > 1 7.2 [ 1− ( − 1 + 3.6min { 2+h 3.6+2h , (h+2)(h2+4) 7.2h−4h2+3.6h3 })2 ] . Finally, from 1− b2 + b1b3 − b23 > 0, we obtain c1(x ∗)2 + c2x ∗ + c3 > 0, i.e., c1 ( 1 + √ 1− 7.2a 3.6 )2 + c2 ( 1 + √ 1− 7.2a 3.6 ) + c3 > 0, where c1, c2, and c3 are expressed as (7). As a simple illustration, consider parameter choice (a, h) = (−2, 0.01). Substituting into the Jury conditions yields 1 + b1 + b2 + b3 ≈ 1.6 × 10−6 > 0, 1 − b1 + b2 − b3 ≈ 7.94 > 0, 1 − b2 + b1b3 − b23 ≈ 1.3 × 10−6 > 0, and |b3| ≈ 0.986 < 1. All conditions are satisfied, confirming that E∗ 1 is stable. In contrast, for (a, h) = (0.01, 0.01) the third M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 6 of 27 inequality becomes negative, which indicates that E∗ 1 is unstable. These numerical results are consistent with the stability regions shown in Fig. 1. In the sense of linearization, the stability conditions derived (6) are local in the sense of linearization: They characterize the behavior in the neighborhood of the equilibrium using the Jacobian and Jury criteria. They do not establish global asymptotic stability nor do they describe transient behavior (e.g., basins of attraction, multistability, or transient chaos). Consequently, the numerical explorations that follow are intended to complement the local theory rather than imply global guarantees. After obtaining the analytical conditions through the Jury Criteria, we numerically verified the results by calculating the Jacobian eigenvalues on the parameter grid (a, h). (a) (b) Figure 1: Two-parameter stability maps for the standard Euler scheme in the (a, h) plane. Panel (a) shows the equilibrium E∗ 1 and panel (b) shows the equilibrium E∗ 2 . Blue indicates stability and gray denotes instability. Regions with a > 1/7.2 (to the right of the dotted vertical line) have no equilibria. Figure 1 shows the numerical stability map for the equilibrium points E∗ 1 and E∗ 2 using the standard Euler scheme. Consistent with the analytical results, only E∗ 1 has a stability domain in the region a < 0 with sufficiently small h, whereas E∗ 2 is always unstable for all pairs (a, h). Furthermore, the larger the value of h, the narrower the range of a that results in stability. 2.2. Vanishing step-size limit (h → 0) The equilibrium of the discrete-time system (2) coincides with that of the continuous- time system (1). Moreover, the Jacobian J(E) can be written as J(E) = I + hA(E), M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 7 of 27 where A(E) is the Jacobian of ODE (1) at equilibrium E = (x, 0, x). This representation allows a direct comparison of the spectral properties of discrete and continuous systems. Theorem 2 (Spectral consistency as h → 0). Let µ1, µ2, µ3 be the eigenvalues of A(E) at an equilibrium E of (1), and let λ1(h), λ2(h), λ3(h) be the eigenvalues of J(E) at the same point. Then λi(h) = 1 + hµi +O(h2), h → 0, i = 1, 2, 3. In particular, if ℜ(µi) < 0 for all i, then there exists h∗ > 0 such that |λi(h)| < 1 for all 0 < h < h∗. Proof. The characteristic equation of J(E) is det(λI − J(E)) = 0, that is, det((λ − 1)I − hA(E)) = 0. Let ν = (λ − 1)/h. Then ν satisfies det(νI − A(E)) = O(h), Thus, ν = µi +O(h) for some eigenvalues µi in A(E). Hence λ = 1+ hµi +O(h2). If ℜ(µi) < 0, then for sufficiently small h one has |1 + hµi| < 1, which indicates the stability of the discrete map. Corollary 1. If an equilibrium E of the ODE (1) satisfies the Routh–Hurwitz conditions (i.e., all eigenvalues µi of A(E) have negative real parts), then the corresponding equilib- rium of the Euler map (2) is locally asymptotically stable for sufficiently small h > 0. Proof. From Theorem 2, the eigenvalues of the Euler Jacobian are λi(h) = 1 + hµi + O(h2) as h → 0. If ℜ(µi) < 0, then |1 + hµi| < 1 for h small enough, and thus each Jury condition (5) is satisfied. Equivalently, one can check directly that the Jury inequalities reduce to the Routh–Hurwitz inequalities for A(E) in the limit h → 0. Thus stability of the continuous-time system implies stability of the discretized system for sufficiently small h. Theorem 2 and Corollary 1 provide several important insights into the relationship between the discrete Euler map and continuous time dynamics. An immediate consequence of Theorem 2 is that the stability of the Euler scheme depends explicitly on step size h. Because the discrete eigenvalues satisfy λ(h) = 1 + hµ, stability requires |1 + hµ| < 1 for each eigenvalue µ of the continuous Jacobian. For a real negative eigenvalue µ = −α with α > 0, this inequality reads |1 − hα| < 1, which simplifies to −1 < 1 − hα < 1 hence, 0 < h < 2/α. For a general complex eigenvalue µ = σ + iω with σ < 0, we compute |1 + hµ|2 = (1 + hσ)2 + (hω)2 = 1 + 2hσ + h2|µ|2. Therefore, the condition |1 + hµ| < 1 is equivalent to 1 + 2hσ + h2|µ|2 < 1, or 2hσ + h2|µ|2 < 0. Because h > 0 and σ < 0, this inequality is reduced to 0 < h < −2σ/|µ|2. In other words, the admissible step size shrinks with both the size of the negative real part, σ and M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 8 of 27 the magnitude of the imaginary part, ω. A conservative guideline that applies uniformly across all eigenvalues is therefore 0 < h < 2 maxi |ℜ(µi)| . This bound highlights the intuitive fact that the stronger the damping in the continuous system (the larger |ℜ(µi)|), the smaller is the admissible step size for the Euler scheme. For example, when a = −2, the dominant eigenvalue has ℜ(µ) ≈ −1.16; therefore, the bound yields h < 1.72, and the choice h = 0.01 used in our simulations lies safely within the stable regime. Beyond the step size restriction, another important implication con- cerns the qualitative changes that arise when h becomes too large. Mathematically, such spurious bifurcations occur when an eigenvalue of the discrete Jacobian J(E) = I+hA(E) crosses the unit circle, that is, when |1 + hµ| = 1 for some eigenvalue µ of the continuous Jacobian. At this point, the equilibrium of the discrete system loses stability even though the continuous system remains stable. For example, if µ = −α < 0, then λ = 1 − hα. When h = 2/α, the discrete eigenvalue becomes λ = −1, producing a change in stability that mimics period-doubling bifurcation. However, at the continuous level, no parameter variation has occurred, Therefore, this bifurcation is purely numerical. More generally, for µ = σ + iω with σ < 0, the crossing condition (1 + hσ)2 + (hω)2 = 1 indicates that a pair of complex-conjugate eigenvalues of the Euler map lies on the unit circle at points different from ±1, which numerically resembles a Neimark–Sacker bifur- cation. Again, the continuous system remained stable and the observed transition was entirely artificial. As h → 0, all such spurious bifurcations vanish, leaving only genuine bifurcations governed by the physical parameter a. Hence, bifurcation diagrams with re- spect to h should be understood not as reflecting true dynamical features of the ODE, but rather as a measure of the sensitivity and limitations of the Euler discretization. 2.3. Numerical Simulations and Dynamical Behavior Having established the analytical stability conditions and their consistency with the continuous-time flow as h → 0, we complement the theoretical results with the numerical simulations. These simulations are designed not only to reveal the global organization of the dynamics through bifurcation diagrams but also to illustrate the qualitative behavior of trajectories by means of representative phase portraits. The initial condition (0, 1.6, 3) was selected because it yields bounded trajectories and representative dynamical behaviors (equilibria, periodic orbits, and chaos) across the chosen parameters, in line with previous studies on quadratic systems without equilibrium [14]. To ensure that only the asymptotic behavior of the system was analyzed, each simulation was run for 40,000 iterations, and the first 35,000 were discarded as transient. This procedure eliminates the influence of the initial settling phase and guarantees that the reported results reflect long-term dynamics (equilibria, periodic orbits, and chaotic attractors). M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 9 of 27 One-parameter bifurcations. System (2) depends on the parameters a and h, with equilibrium E = (x, 0, x). First, we examine the one-parameter bifurcations using the local maxima xmax of xt. Bifurcation diagrams were generated by scanning the con- trol parameters over a uniform grid and plotting the numerically detected maxima. The accompanying largest Lyapunov exponent was computed using the standard QR-based algorithm to confirm the transition between periodic and chaotic regimes. As shown in Figure 2, the bifurcation diagram with respect to a exhibits a typical period doubling route to chaos. Figure 3 illustrates the effect of step size h, where increasing h eventually destabilizes the dynamics and induces numerical chaos in agreement with the theoretical step-size restriction. (a) (b) Figure 2: One-parameter bifurcation (standard Euler) with respect to a at fixed h = 0.01. (a) Local maximum xmax of xt. (b) Lyapunov spectrum (λ1, λ2, λ3). Each run used 40,000 iterations, with the first 35,000 discarded as transients. Positive values of the largest Lyapunov exponent indicate chaotic intervals, while clustered xmax bands (panel a) correspond to periodic windows. The largest Lyapunov exponent was computed following the standard algorithm of Wolf et al. [34], in which the separation of nearby trajectories was measured iteratively and renormalized along the orbit to avoid numerical overflow. Specifically, for each parameter set, we evolved the system for 5× 104 iterations, discarded the first 104 as transients, and applied a renormalization procedure at each step to estimate the average exponential rate of divergence. A positive value of the computed exponent is considered as an indicator of chaos. This method is widely used in the study of discrete dynamical systems and ensures the consistency of chaotic diagnostics (see also [35, 36]). M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 10 of 27 (a) (b) Figure 3: One-parameter bifurcation (standard Euler) with respect to step size h at fixed a = 1. (a) Local maxima xmax; (b) Lyapunov spectrum (λ1, λ2, λ3). Each run uses 40,000 iterations with the first 35,000 discarded as transients. The spectrum highlights chaotic regimes (positive λ1) and step-size–induced numerical instabilities. Two-parameter bifurcation. To explore the combined effects of a and h, we con- structed a two-parameter diagram in the (a, h) plane. The maximum value xmax in the trajectory was recorded, and the largest Lyapunov exponent was computed in parallel. The resulting diagrams in Figure 4 provide a compact global picture of the dynamics, including the periodic windows, embedded chaotic regions, and divergence boundaries. Attractor dimension Beyond the largest Lyapunov exponent, we also quantified at- tractor complexity by computing the Kaplan–Yorke (KY) dimension from the full Lya- punov spectrum. For a typical 3D chaotic map, we find λ1 > 0, λ2 ≈ 0, λ3 < 0, and thus j = 2 and DKY = 2 + (λ1 + λ2)/|λ3|. When trajectories converge to an equilibrium or periodic orbit, the spectrum collapses and DKY takes integer values (0 or 1), whereas in chaotic regimes it lies strictly between 2 and 3. The variation of DKY with respect to a and h is shown in Figure 5, confirming the fractal geometry of chaotic attractors. M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 11 of 27 (a) (b) Figure 4: Two-parameter bifurcation diagrams in the (a, h) plane for the standard Euler scheme. (a) Map of xmax: colour bar encodes the magnitude of the largest x values; black regions correspond to divergence. (b) Largest Lyapunov exponent map computed via QR decomposition: negative values (yellow–magenta shades) correspond to stable or periodic dynamics, while positive values (green–blue shades) indicate chaotic regimes. (a) (a) (b) (b) Figure 5: Kaplan–Yorke dimension DKY for the Euler discretization. (a) Variation in DKY with parameter a at a fixed h = 0.01. (b) Variation in DKY with step size h at a fixed a = 1. Stable equilibria and periodic orbits yield integer DKY values, while chaotic regimes are characterized by noninteger dimensions typically between 1 and 3. Phase portraits Guided by bifurcation analysis, we illustrate the representative tra- jectories. The phase portraits were generated by iterating the system (2) with step size h = 0.01 for 60,000 steps, starting from the initial condition (x0, y0, z0) = (0, 1.6, 3). Only M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 12 of 27 the last 5,000 points were plotted to remove the transients. For a = −2, the system con- verges to a stable equilibrium, whereas for a = 1 it exhibits a strange attractor, consistent with previous studies [14, 37–39]. Representative trajectories are illustrated in Figure 6, where both stable and chaotic behaviors are clearly distinguished. The chaotic attractor displays a folded-band structure characteristic of quasiperiodic breakdown. (a) (b) (c) (d) Figure 6: Phase portraits of system (2) with h = 0.01 and initial condition (0, 1.6, 3). Panels (a)–(b): parameter a = −2, the trajectory converging to the stable equilibrium E∗ 1 (red point), and E∗ 2 (blue point) remain unstable. Panels (c)–(d): Parameter a = 1, showing a chaotic attractor with a folded-band geometry. Each plot was based on 40,000 iterations, and the first 35,000 transients were discarded. Panels (a) and (c) display the three-dimensional trajectories, whereas panels (b) and (d) show the (x, y) projections. M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 13 of 27 3. Discrete-Time System via Nonstandard Mickens Scheme We follow one of the nonstandard difference schemes proposed by Mickens [32] where the term on the right-hand side, which consists of the same variable on the left-hand side with a negative coefficient, is replaced by −1.8xtzt = −1.8xtzt+1. This modification follows Mickens’ principle of nonlocal approximations, in which nonlinear terms are discretized asymmetrically to better preserve the dynamics of the continuous system. Therefore, we have a new system as follows xt+1 = xt + hyt yt+1 = yt + h(−xt + zt) zt+1 = [zt + h(zt − 2xtyt − a)]/[1 + 1.8hxt]. (9) 3.1. Equilibrium and Stability Criteria It is easy to verify that the nonstandard system (9) has the same equilibrium points as the system (2). The Jacobian matrix of nonstandard system (9) evaluated at any point E = (x, 0, x) ∈ R3 is JNS(E) =  1 h 0 −h 1 h 1.8h((a−x)h−x) (1.8hx+1)2 − 2hx 1.8hx+1 1+h 1.8hx+1  (10) The characteristic polynomial of (10) is pNS(λ) = λ3 + b1λ 2 + b2λ+ b3 where b1 = −3.6hx+ h+ 3 1.8hx+ 1 b2 = 1.8h3x+ 2hx2 + h2 + 1.8hx+ 2h+ 3 1.8hx+ 1 b3 = −1− 1.8a h4 + ( −3.6x2 − 1 ) h3 + (−1− 3.8x)h2 + (−1− 1.8x)h (1 + 1.8hx)2 Compared with the standard Euler scheme, the coefficients of the characteristic polynomial are more complicated owing to the rational dependence on both h and x introduced by the nonstandard denominator. Based on the first Jury’s condition in (5), we have 1 + b1 + b2 + b3 = ( 16.2hx2 + 18x+ 9a ) h3 (1 + 1.8hx)2 > 0 ⇒ ( 16.2hx2 + 18x+ 9a ) > 0, (11) M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 14 of 27 Note that equation (11) is quadratic concerning x such that it is positive definite, the coefficient of x2 must be positive, and the discriminant must be negative. This is equivalent to a < − 1 1.8h . (12) Next, for the second Jury’s condition, we obtain the following equation 1− b1 + b2 − b3 = ( 20.16h2 + 3.24h4 ) x2 + ( 25.6h+ 3.6h3 ) x+ 8 + 2h2 + 1.8h3a (1.8hx+ 1)2 > 0 ⇒ ( 20.16h2 + 3.24h4 ) x2 + ( 25.6h+ 3.6h3 ) x+ 8 + 2h2 + 1.8h3a > 0. (13) As in the case of the first condition, it is sufficient to consider the coefficient of x2 positive and the discriminant negative. This is equivalent with a > −405h4 − 2520h2 + 320 729h5 + 4536h3 . (14) For the third Jury’s condition 1−b2+b1b3−b3 2 > 0, this inequality cannot be simplified in closed form, but it can be checked numerically for the given parameter values. c1x 4 + c2x 3 + c3x 2 + c4x+ c5 > 0 (15) where x = x∗ c1 = −31.4928h6 − 68.9472h4 c2 = −75.816h5 − 179.64h3 c3 = −12.96h5a− 75.24h4 − 174.64h2 c4 = −13.68ah4 − 34.6h3 − 75h c5 = −12− 3.24h6a2 − 3.6h5a− h4 − 3.6h3a− 6h2. Finally, for the last Jury’s condition b2 > 3, we have 1.8h3x+ h2 − 1.6hx 1.8hx+ 1 > 0 (16) Note that the numerator of the equation (16) is positive definite for h > 0 so the sign of the equation (16) is affected by the denominator only, namely when the condition x∗ > − 1 1.8h . (17) Note that equilibrium E∗ 1 is always satisfied (17) where for E∗ 2 the (17) is equivalent with a < 1 7.2 − 1 7.2 ( 1 h + 1 )2 . (18) Therefore, the two equilibria of the non-standard system (9) are asymptotically stable if the conditions in (12), (14), (15), and (17) are satisfied. M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 15 of 27 Theorem 3. Consider the nonstandard system (9) with equilibrium points E∗ 1,2 = ( 1± √ 1− 7.2a 3.6 , 0, 1± √ 1− 7.2a 3.6 ) . The equilibrium points E∗ 1,2 are locally asymptotically stable if the following conditions hold.  −405h4 − 2520h2 + 320 729h5 + 4536h3 < a < − 1 1.8h , x∗ > − 1 1.8h , c1(x ∗)4 + c2(x ∗)3 + c3(x ∗)2 + c4x ∗ + c5 > 0, (19) where c1 = −31.4928h6 − 68.9472h4, c2 = −75.816h5 − 179.64h3, c3 = −12.96h5a− 75.24h4 − 174.64h2, c4 = −13.68ah4 − 34.6h3 − 75h, c5 = −12− 3.24h6a2 − 3.6h5a− h4 − 3.6h3a− 6h2. As in the Euler case (Section 2), these stability results are purely local; global stability and transient dynamics are outside the present scope. Furthermore, using the same pro- cedure, we also evaluate the stability of the equilibrium points E∗ 1 and E∗ 2 in the Mickens nonstandard scheme. The two-parameter stability map in the (a, h) plane is shown in Figure 7, which is constructed based on the spectral radius of the non-standard Jacobian. Figure 7 shows the numerical results for the non-standard scheme. The equilibrium point E∗ 1 remains stable at the interval a < 0; however, its stability region is more limited than that of the standard scheme because of the additional condition 1 + 1.8hx∗ > 0. Meanwhile, although the Jury conditions indicate that both E∗ 1 and E∗ 2 could in principle be stable under suitable (a, h) pairs, the numerical exploration shows that E∗ 2 possesses no significant stability domain within the investigated ranges. In practice, only E∗ 1 exhibits a significant stability region. This highlights that the main difference between the standard and nonstandard schemes emerges more clearly when h increases. 3.2. Vanishing step-size limit (h → 0): Nonstandard Mickens In the nonstandard scheme (9), the expansion of the denominator (1 + 1.8hx)−1 = 1− 1.8hx+O(h2) (1 + 1.8hx)−2 = 1− 3.6hx+O(h2) shows that the Jacobian matrix at the equilibrium point E = (x∗, 0, x∗) can be written in the form JNS(E) = I + hA(E) +O(h2), h → 0, M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 16 of 27 (a) (b) Figure 7: Two-parameter stability maps for the nonstandard Mickens scheme in the (a, h) plane. Panel (a): Equilibrium E∗ 1 ; panel (b): equilibrium E∗ 2 . Blue regions denote stable parameter pairs, gray regions indicate unstable equilibria, and white regions correspond to parameter values in which equilibrium does not exist. Stability was determined from the spectral radius of JNS(E) at equation (10). where A(E) is the Jacobian of the ODE (1) at the equilibrium E = (x, 0, x). This confirms that the spectral behavior of the nonstandard scheme is consistent with that of the ODE system. Theorem 4 (Spectral consistency, nonstandard). If µ1, µ2, µ3 are eigenvalues of A(E), while λi(h) are eigenvalues of JNS(E), then λi(h) = 1 + hµi +O(h2), h → 0, i = 1, 2, 3. A direct consequence of this consistency is the connection between the Routh–Hurwitz condition for continuous systems and Jury stability for nonstandard maps. Corollary 2. If all ℜ(µi) < 0, then there exists h∗ > 0 such that |λi(h)| < 1 for every 0 < h < h∗. In other words, local asymptotic stability in the continuous system guarantees the stability of the equilibrium point in the non-standard scheme when h is sufficiently small. Thus, for h approaching zero, the nonstandard scheme not only shares equilibrium points with the continuous system but also has consistent local stability properties. This clarifies that the main difference with the standard scheme arises only for a larger h, and not for small limits. Remark 1 (Domain and step size conditions). In addition to the spectral conditions above, the non-standard scheme requires 1+ 1.8hx∗ > 0 for the denominator to be non-zero. this condition is satisfied automatically for h → 0. Furthermore, the approximation |1+hµ| < 1 M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 17 of 27 provides a conservative bound for h: if µ = −α < 0 then 0 < h < 2/α, while if µ = σ+ iω with σ < 0 then 0 < h < −2σ |µ|2 . As a uniform guideline, the following applies 0 < h < 2 maxi |ℜ(µi)| , which explains the reason for choosing a small h in the numerical simulations, while also warning about the possibility of pseudo-bifurcation if h is too large. 3.3. Numerical Exploration and Comparison with the Standard Scheme Having established that the nonstandard scheme is consistent with the continuous sys- tem as h → 0, We now present numerical explorations of the same order as the standard scheme. one-parameter bifurcations, two-parameter diagrams, attractor dimension analy- sis, representative phase portraits, and a final comparison with the continuous model. One-Parameter Bifurcation. Figure 8 shows the bifurcation profile and corresponding largest Lyapunov exponent as the parameter a varies with fixed step size h = 0.01. Both panels are consistent and exhibit the same transition structure. Similarly, Figure 9 presents the one-parameter diagrams with varying step size h and fixed a = 1. The xmax map and the Lyapunov exponent corroborate each other, illustrating how decreasing h affects the onset of complex dynamics. This behavior parallels the Euler scheme, although the nonstandard denominator delays the onset of spurious bifurcations. Two-parameter diagrams in (a, h). To explore the joint influence of a and h, Fig- ure 10 shows the two-parameter diagram based on xmax together with the corresponding largest Lyapunov exponent. The peak-based map organizes periodic windows and thin transition bands, whereas the Lyapunov map separates the stable (negative) and chaotic (near-zero/positive) regions. Compared with the Euler case, the stability region of E∗ 1 is narrower and E∗ 2 does not exhibit a significant stability window. M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 18 of 27 (a) (b) Figure 8: One-parameter bifurcation of the nonstandard Mickens scheme with respect to a at h = 0.01. (a) Local maxima xmax after discarding the transients. (b) Lyapunov spectrum (λ1, λ2, λ3), confirming regular-to-chaotic transitions, where positive values of λ1 indicate chaos, while negative values correspond to periodic or stable regimes. Each simulation used 40,000 iterations with the first 35,000 discarded as transients. (a) (b) Figure 10: Two-parameter diagrams in the (a, h) plane for the nonstandard Mickens scheme. (a) Map of xmax: the color bar encodes the magnitude of local maxima, where different color bands correspond to different orbit periods, and black regions denote di- vergence. The white regions indicate either convergence to a fixed point (no oscillatory maxima recorded) or parameter values where no equilibrium exists (for a > 1/7.2). (b) The Largest Lyapunov exponent map computed via QR decomposition: negative val- ues (magenta shades) correspond to periodic or stable dynamics, whereas positive values (green–blue shades) indicate chaotic regimes. Each grid point was iterated for 40,000 steps with the first 35,000 discarded as transients. M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 19 of 27 (a) (b) Figure 9: One-parameter bifurcation of the nonstandard Mickens scheme with respect to h at a = 1. (a) Local maxima xmax after discarding the transients. (b) Lyapunov spectrum (λ1, λ2, λ3), where positive values of λ1 indicate chaos, and negative values correspond to periodic or stable dynamics. Each simulation used 40,000 iterations with the first 35,000 discarded as transients. Attractor dimension. The Kaplan–Yorke (KY) dimension was computed for the Mick- ens scheme, using the same procedure described in the Euler case. Figure 11 shows that DKY = 0 when the trajectories converge to equilibrium and integer values (1, 2) for peri- odic orbits, as also seen in the bifurcation diagrams. Occasionally, intermediate non-integer values between 1 and 2 appear, corresponding to quasiperiodic or transition states. Fully chaotic regimes are characterized by 2 < DKY < 3, which is consistent with the positive Lyapunov exponents and broadband bifurcation structures. The results closely mirror those of the Euler scheme, with both giving DKY = 0 in stable equilibrium, integer values for periodic dynamics, 1 < DKY < 2 in quasiperiodic regimes, and 2 < DKY < 3 in chaotic regimes. Representative phase portraits. Figure 12 shows the typical phase portraits for h = 0.01 and the initial condition (0, 1.6, 3). For a = −2 trajectories converge to a stable equilibrium, whereas for a = 1 a chaotic attractor emerges, closely mirroring the standard scheme with minor differences in the chaotic regime. The attractor geometry appears slightly sharper than in Euler due to the rational structure of the discretization. Comparison with the continuous system. Finally, Figure 13 compares the xmax- based bifurcation curves of the continuous system (black, ode45) with the standard Euler (red) and Mickens non-standard (blue) discretizations at h = 0.01. For a < 0 the three agree closely; near the transition around a ≈ 1 both discretizations capture the onset of M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 20 of 27 (a) DKY vs. a at fixed h = 0.01 (b) DKY vs. h at fixed a = 1 Figure 11: Kaplan–Yorke dimension DKY for the nonstandard Mickens discretization. (a) Variation of DKY with parameter a at fixed h = 0.01. (b) Variation in DKY with step size h at a fixed a = 1. Stable equilibria give DKY ≈ 0, periodic orbits yield integer values (e.g., DKY = 1, 2), and chaotic regimes are characterized by non-integer values 2 < DKY < 3. Each computation was performed with 40,000 iterations after discarding the first 35,000 transients. complex dynamics with small amplitude biases relative to the continuous case; for a > 1 all three show a chaotic spread, with the non-standard scheme, producing broader chaotic bands while maintaining the qualitative structure. In addition to the comparison with the continuous system, it is useful to highlight the direct differences between the Euler and Mickens discretizations, which are reported in Table 1. Both schemes reproduce the qualitative features of the ODE for small h, but the Mickens scheme shows greater robustness in preserving attractor dimensions and delays numerical instabilities at larger h. Furthermore, it is worth situating the NSFD approach within a broader landscape of numerical methods. Classical high-order schemes, such as Runge–Kutta (e.g., ode45) and semi-implicit integrators, are widely used for their accuracy when the step size is small. However, these conventional methods do not necessarily preserve qualitative features, such as positivity, boundedness, or equilibrium stability, when larger step sizes are employed. In contrast, Mickens’ nonstandard scheme was explicitly constructed to maintain such structural properties, which explains its robustness in capturing chaotic and bifurcating dynamics under coarse discretizations. This broader perspective emphasizes that NSFD schemes provide a valuable complementary alternative to classical integrators for studying discrete nonlinear dynamics. M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 21 of 27 (a) (b) (c) (d) Figure 12: Phase portraits of the nonstandard system (9) with h = 0.01 and initial condition (0, 1.6, 3). Panels (a)–(b): parameter a = −2, the trajectory converging to the stable equilibrium E∗ 1 (red point), and E∗ 2 (blue point) remain unstable. Panels (c)– (d): Parameter a = 1, showing a chaotic attractor with a folded-band geometry. Each plot is based on 40,000 iterations, discarding the first 35,000 as transients. Panels (a) and (c) display three-dimensional trajectories, whereas panels (b) and (d) show the (x, y) projections. 4. Conclusion This study compared two discrete-time versions of a three-dimensional nonlinear sys- tem: one discretized by the standard forward Euler method and the other by a nonstandard Mickens scheme. Both approaches share the same equilibrium points, but their stability M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 22 of 27 Figure 13: Comparison of xmax-based bifurcation diagrams of the continuous system (black, computed with ode45), standard Euler discretization (red), and Mickens non- standard scheme (blue), all with a step size h = 0.01. For small h, both discretizations reproduce the continuous dynamics, but near the onset of chaos (a ≈ 1), the nonstandard scheme remains closer to the continuous reference while the Euler method shows larger deviations. properties differ. We derived discrete systems, analyzed their equilibria, and applied the Jury criterion to determine local stability. The results show that the Euler scheme sta- bilizes only one equilibrium. In the nonstandard scheme, the theoretical conditions allow the possibility of stabilizing both equilibria; however, numerical exploration indicates that only E∗ 1 exhibits a robust stability domain, whereas E∗ 2 remains practically unstable across most parameter ranges. For the Euler discretization, we derived the Jacobian characteristic polynomial and applied the Jury criterion. The results show that equilibrium E∗ 1 is stable for certain (a, h) values, whereas E∗ 2 is always unstable. The stability region, expressed through inequalities in the parameters, was numerically confirmed. The nonstandard Mickens scheme modifies the nonlinear term to preserve the system structure. Jacobian and characteristic polynomials were derived, and the Jury criterion yielded nonlinear conditions in (a, h). Unlike Euler, both equilibria E∗ 1 and E∗ 2 can be stable under appropriate parameters. A set of numerical simulations were performed to delve deeper into the dynamic conse- quences of these theoretical findings. Phase portraits illustrated that both methods could replicate the periodic and chaotic behaviors observed in the continuous system, although differences appeared in the detailed structure of the chaotic attractors. The one-parameter bifurcation diagrams of a and h show that both discretizations display complex nonlinear characteristics including bifurcations, periodic windows, and chaotic zones. Nonetheless, M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 23 of 27 Table 1: Comparison between Euler and Mickens discretizations. Aspect Euler scheme Mickens scheme Equilibria Shares with continuous sys- tem Shares with continuous sys- tem Stability region Larger for E∗ 1 and E∗ 2 at small h Narrower; E∗ 2 lacks significant stable window Bifurcation onset Spurious transitions appear earlier as h grows Onset of instability delayed by nonstandard denominator Lyapunov exponents Detect chaos reliably for small h Same reliability; slightly smoother near thresholds Kaplan–Yorke dimen- sion 2 < DKY < 3 in chaos, col- lapses to 2 in stability Same, with better robustness at larger h Phase portraits Qualitatively correct; slightly blurred chaotic attractor Qualitatively correct; attrac- tor geometry sharper Agreement with ODE Good at small h, deteriorates faster for large h Good at small h, remains closer for moderate h the nonstandard scheme demonstrated superior numerical stability, particularly for larger h values, where the standard Euler method began to deviate significantly or failed to maintain bounded trajectories. Two-parameter bifurcation diagrams and Lyapunov exponent landscapes in the (a, h) plane offer a more profound understanding of the overall behavior of each discrete model. In the standard scheme, chaos predominantly appears in areas where a > −1 and h > 0.2, as suggested by the Lyapunov exponents that are near zero or positive. Conversely, the nonstandard scheme displayed a more complex pattern: chaotic regions emerged even with relatively low h values, whereas stable areas were maintained more extensively across the parameter space. Significantly, the Lyapunov exponents from the nonstandard scheme span a broader range, underscoring the existence of more complex dynamics and shifts between stability and chaos. These findings highlight the benefits of employing a nonstandard approach for dis- cretizing nonlinear systems with sensitive dynamics. The Mickens-based technique not only maintains the equilibrium structure but also enhances stability and preserves chaotic dynamics over a wider range of parameter values and time step sizes. This makes it an excellent option for long-term, coarse-grained numerical experiments. These results also emphasize the significance of choosing suitable numerical schemes when studying discrete- time dynamical systems, particularly for capturing intricate structures such as bifurcations and chaotic attractors. Future studies might expand these methods to include delayed, fractional, or stochastic versions of the model, or explore the effects of other nonstandard discretization rules on various types of dynamical systems. Although the findings are encouraging, this study had several limitations. First, the stability analysis was performed by linearizing around fixed points, which only describes local dynamics and does not consider global phenomena, such as basin structures, co- M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 24 of 27 existing attractors, multi-stability, invariant manifolds, or transient chaos. Therefore, a comprehensive global study is necessary to fully characterize discrete dynamics. Second, the parameter space examined, although broad, was restricted to two control parameters (a and h); other factors such as perturbations, delay effects, or coupling with additional vari- ables were not included. Third, while the Mickens scheme used here offers improvements over the Euler method, other types of nonstandard finite-difference schemes with different denominator functions or implicit forms may provide even more precise long-term behav- ior, especially in stiff systems. Finally, the present work remains mathematical in nature; embedding the model into physical or biological contexts, for example, chaotic electronic circuits, fractional-order epidemiological transmission models, or ecological predator–prey systems—could provide both practical applications and validation of the proposed dis- cretization strategies. Future studies should extend these results in several ways. Possible paths include higher-order nonstandard or symplectic schemes that preserve geometric invariants; appli- cations to systems with delays, fractional derivatives, or external forcing; and global bi- furcation analyses, including invariant manifolds and basin boundaries. Computationally, adaptive step-size methods or machine learning-assisted bifurcation tracking may improve efficiency. Finally, stochastic extensions of discretized systems can clarify the role of noise in discrete chaos. In particular, it would be fruitful to apply the present discretization strategies to concrete models such as Lorenz-like chaotic systems, fractional-order epi- demiological models, or delayed feedback systems, which would provide both theoretical challenges and practical relevance for future exploration. We hope that these directions will stimulate further studies in discrete nonlinear dynamics. Acknowledgements This research was funded by the Faculty of Science and Mathematics, Universitas Diponegoro, with funding from Selain APBN Universitas Diponegoro 2025 (Grant Number 26.G/UN7.F8/PP/II/2025). References [1] M. Molaie, S. Jafari, J. C. Sprott, and S. M. R. H. Golpayegani. Simple chaotic flows with one stable equilibrium. International Journal of Bifurcation and Chaos, 23(11), 2013. [2] V. T. Pham, S. Jafari, C. K. Volos, and S. Vaidyanathan. Hidden attractors in a chaotic system with an exponential nonlinear term. The European Physical Journal Special Topics, 224(8):1507–1517, 2015. [3] F. R. Tahir, S. Jafari, V. T. Pham, C. Volos, and X. Wang. A novel no-equilibrium chaotic system with multiwing butterfly attractors. International Journal of Bifur- cation and Chaos, 25(4), 2015. [4] D. Dudkowski, S. Jafari, T. Kapitaniak, N. V. Kuznetsov, G. A. Leonov, and M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 25 of 27 A. Prasad. Hidden attractors in dynamical systems. Physics Reports, 637:1–50, 2016. [5] M. Parsamanesh and M. Gümüş. Qualitative study for the system of waste plastic management in the ocean: A discrete-time deterministic model. Communications in Nonlinear Science and Numerical Simulation, page 108617, 2025. [6] D. T. Wood, H. V. Kojouharov, and D. T. Dimitrov. Universal approaches to approxi- mate biological systems with nonstandard finite difference methods. Mathematics and Computers in Simulation, 133:337–350, 2017. [7] M. Gümüş and K. Türk. A note on the dynamics of a covid-19 epidemic model with saturated incidence rate. The European Physical Journal Special Topics, pages 1–8, 2024. [8] J. Nainggolan, M. F. Ansori, and H. Tasman. An optimal control model with sen- sitivity analysis for covid-19 transmission using logistic recruitment rate. Healthcare Analytics, page 100375, 2025. [9] Y. Halim, A. Allam, M. Fečkan, I. Redjam, and M. Gümüş. Exploring competitive, mutualistic, and other interactions in a discrete community model. Journal of Applied Mathematics and Computing, pages 1–30, 2025. [10] M. F. Ansori, S. Brianzoni, and G. Campisi. Bifurcations and complex dynamics in a banking duopoly model with macroprudential policy. Physica A: Statistical Mechanics and its Applications, 641:129730, 2024. [11] M. F. Ansori, N. Y. Ashar, and H. K. Fata. Logistic map-based banking loan dynamics with central bank policies. Journal of Applied Nonlinear Dynamics, 14(03):561–574, 2025. [12] S. Lasfar, M. Elkarmouchi, K. Hattaf, and N. Yousfi. Mathematical analysis of the dynamics of a fractional economic cycle model and the existence of solutions by means of fixed point theory. Advances in Fixed Point Theory, 14, 2024. [13] I. Ispolatov, S. Allende, M. Doebeli, and V. Madhok. Chaos in high-dimensional dissipative dynamical systems. Scientific Reports, 5(1):12506, 2015. [14] S. Jafari, J. C. Sprott, and S. M. R. H. Golpayegani. Elementary quadratic chaotic flows with no equilibria. Physics Letters A, 377(9):699–702, 2013. [15] A. Akgul and I. Pehlivan. A new three-dimensional chaotic system without equi- librium points, its dynamical analyses and electronic circuit application. Technical Gazette, 23(1):209–214, 2016. [16] J. O. Maaita, Ch. K. Volos, I. M. Kyprianidis, and I. N. Stouboulos. The dynamics of a cubic nonlinear system with no equilibrium point. Journal of Nonlinear Dynamics, 2015:257923, 2015. [17] Y. Gong, X. Wang, and Y. Liu. New 4d chaotic system with hidden attractors and self-excited attractors and its application in image encryption based on rng. Physica A: Statistical Mechanics and its Applications, 601:127–145, 2022. [18] P. Djorwe, E. Ngouonkadi, and P. Woafo. Hidden attractors and metamorphoses of basin boundaries in optomechanics. Nonlinear Dynamics, 111:1555–1572, 2023. [19] N. Kuznetsov, G. Leonov, V. Bragin, et al. Hidden attractors in chua circuit: math- ematical theory meets physical experiments. Nonlinear Dynamics, 111:2575–2591, M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 26 of 27 2023. [20] X. Yang and C. Lai. Construction and implementation of discrete memristive hyperchaotic map with hidden attractors and self-excited attractors. Integration, 89:155–166, 2023. [21] J. Maaita and D. Prousalis. A comparison between four chaotic indicators in sys- tems with hidden attractors. Journal of Computational and Nonlinear Dynamics, 20(1):011008, 2024. [22] H. Zhang, J. Li, and G. Chen. Self-excited and hidden multi-scroll attractors in a novel extremely chaotic system and its cryptographic applications. Physica A: Statistical Mechanics and its Applications, 643:129875, 2025. [23] A. Ouannas, V. T. Pham, A. A. Khennaoui, G. Grassi, and S. Momani. Chaos and control of a three-dimensional fractional order discrete-time system with no equilib- rium and its synchronization. American Institute of Physics Advances, 10(4):045310, 2020. [24] X. Wang, C. Volos, E. Tlelo-Cuautle, V. T. Pham, J. M. Munoz-Pacheco, and S. Ja- fari. A new chaotic system with stable equilibrium: From theoretical model to circuit implementation. Institute of Electrical and Electronics Engineers Access, 5:8851–8858, 2017. [25] J. P. Singh, K. Rajagopal, and B. K. Roy. A new 5d hyperchaotic system with stable equilibrium point, transient chaotic behaviour and its fractional-order form. Pramana, 91(3), 2018. [26] K. Nonlaopon, M. M. Khalsaraei, A. Shokri, and M. Molayi. Approximate solu- tions for a class of predator–prey systems with nonstandard finite difference schemes. Symmetry, 14(8):1660, 2022. [27] R. U. Din, K. A. Khan, A. Aloqaily, N. Mlaiki, and H. Alrabaiah. Using nonstandard finite difference scheme to study classical and fractional order seivr model. Fractal and Fractional, 7(7):552, 2023. [28] K. Kaur and G. Singh. An efficient nonstandard numerical scheme coupled with a compact finite difference method to solve the one-dimensional burgers’ equation. Axioms, 12(6):593, 2023. [29] I. U. Khan, A. Hussain, S. Li, and A. Shokri. Modeling the transmission dynamics of coronavirus using nonstandard finite difference scheme. Fractal and Fractional, 7(6):451, 2023. [30] M. Gümüş and K. Türk. Dynamical behavior of a hepatitis b epidemic model and its nonstandard finite difference scheme. Journal of Applied Mathematics and Comput- ing, 70(4):3767–3788, 2024. [31] M. Hajipour, A. Jajarmi, and D. Baleanu. An efficient nonstandard finite differ- ence scheme for a class of fractional chaotic systems. Journal of Computational and Nonlinear Dynamics, 13(2):021013, 2018. [32] R. E. Mickens. Nonstandard Finite Difference Models of Differential Equations. World Scientific, Singapore, 1994. [33] J. Andaluz, A. A. Elsadany, and G. Jarne. Nonlinear cournot and bertrand-type dy- namic triopoly with differentiated products and heterogeneous expectations. Mathe- M. F. Ansori et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6506 27 of 27 matics and Computers in Simulation, 132:86–99, 2017. [34] A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano. Determining lyapunov exponents from a time series. Physica D: Nonlinear Phenomena, 16(3):285–317, 1985. [35] G. Benettin, L. Galgani, A. Giorgilli, and J. Strelcyn. Lyapunov characteristic ex- ponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. part 1: Theory. Meccanica, 15(1):9–20, 1980. [36] I. Shimada and T. Nagashima. A numerical approach to ergodic problem of dissipative dynamical systems. Progress of Theoretical Physics, 61(6):1605–1616, 1979. [37] H. Broer, C. Simó, and R. Vitolo. Bifurcations and strange attractors in the lorenz-84 climate model with seasonal forcing. Nonlinearity, 15(4):1205–1267, 2002. [38] H. Broer, R. Vitolo, and C. Simó. Quasi-periodic hénon-like attractors in the lorenz-84 climate model with seasonal forcing. In European Conference on Qualitative Theory of Differential Equations 2003, pages 601–606. World Scientific, Singapore, 2005. [39] R. Vitolo, H. Broer, and C. Simó. Quasi-periodic bifurcations of invariant circles in low-dimensional dissipative dynamical systems. Regular and Chaotic Dynamics, 16(2–3):154–184, 2011.