EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
2025, Vol. 18, Issue 3, Article Number 6381
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Analytical and Numerical Investigation of a
Fractional-Order 4D Chaotic System via Caputo

Fractional Derivative

Ilhem Kadri1,∗, Rania Saadeh2, Dalal M. AlMutairi3, Mohammed E. Dafaalla4,
Mohammed Berir5, Mohamed A. Abdoon5

1 Department of Mathematics, University of Oran 1 Ahmed Ben Bella, Oran 31000, Algeria
2 Department of Mathematics, Zarqa University, Zarqa 13110, Jordan
3 Department of Mathematics, College of Science and Humanities, Shaqra University,
Al-Dawadmi 17472, Saudi Arabia
4 Department of Mathematics, College of Science, Qassim University, Buraidah 51452,
Saudi Arabia
5 Department of Mathematics, Faculty of Science, Bakht Al-Ruda University, Duwaym
28812, Sudan

Abstract. This work applies two efficient techniques to investigate fractional-order systems: the
Residual Power Series Method (RPSM) and the Caputo fractional derivative (CFD). Both are
utilized, especially to evaluate chaotic behavior and investigate the complex dynamics associated
with a four-dimensional fractional-order chaotic system. Reliable and efficient numerical simula-
tions with the complexity of chaotic behavior are obtained using the CFD approach. Furthermore,
derived analytical solutions of the fractional-order 4D system using the RPSM are obtained. This
approach is highly preferred as it can handle several beginning conditions, is computationally
efficient, and is numerically stable, hence producing rather precise results. Combining RPSM’s
analytical capacity with CFD’s strong numerical accuracy offers a comprehensive understanding of
the system dynamics. Although RPSM is best for stability and simplicity of use, the CFD technique
is especially helpful because of its great accuracy in predicting chaotic behavior. The suggested
methods precisely define their dynamics, provide exact solutions, and clearly identify chaotic at-
tractors. For instance, representative parameter values used in our analysis are (ν = 0.95, 0.99, 1).
These results show how well CFD- and RPSM-based methods represent and solve challenging
problems in engineering and scientific research.

2020 Mathematics Subject Classifications: 26A33,34H10, 37M05

Key Words and Phrases: Fractional Derivatives, Caputo fractional derivative, Residual Power
Series Method, Chaos, Simulation

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6381

Email addresses: kadri.ilhem@univ-oran1.dz (I. Kadri), rsaadeh@zu.edu.jo (R. Saadeh),
dalmteri@su.edu.sa (D. M. AlMutairi), m.dafaalla@qu.edu.sa (M. E. Dafaalla),
midriss@bu.edu.sa (M. Berir), mabdoon.c@ksu.edu.sa (M. A. Abdoon),

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 2 of 19

1. Introduction

In recent decades, fractional calculus has been an excellent generalization of traditional
integer-order differentiation and integration and has been shown to be a more realistic
description of advanced dynamic phenomena than other traditional descriptions [1–6]. It
is an extremely flexible and realistic mathematical tool for the representation of many
physical, biological, and engineering phenomena [7–10]. Its performance is stronger in
long-term memory and dependency systems, where traditional processes fail [11–14].

Various definitions of fractional derivatives exist, such as those proposed by Riemann-
Liouville, Caputo, Caputo-Fabrizio, and Atangana-Baleanu, which offer different benefits
regarding mathematical precision, ease of calculation, and applicability in general to many
physical systems [15–18]. The clever command of the aforementioned formulations is the
secret to their use in the fields of physics, engineering, and biology, where the incorporation
of fractional calculus is the underlying addition to precision and analytical competence [19–
21]. As classical differentiation finds it difficult to deal with the irregular and infinitely
detailed characters of fractals, fractal derivatives provide a better mathematical instrument
to explain their complicated structures and dynamic behaviors.

The power of fractional-order systems to integrate long-term memory and spatial het-
erogeneity into mathematical models makes them a crucial tool for studying and managing
chaos. This is especially crucial in fields like economics, biology, engineering, and physics,
where system behavior frequently demonstrates historical dependence and fractal features.
Fractional derivatives allow for richer dynamical behavior, such as delayed reactions, on-
going correlations, and more seamless transitions between periodic and chaotic regimes,
by permitting the order of the system to be either constant or variable. Therefore, a key
component in improving our comprehension of complicated systems is fractional calculus.

The Caputo derivative offers a realistic representation of physical systems with memory
effects and is therefore particularly well-adapted to represent real-world dynamics. Its
non-singular kernel can simulate processes with exponential decay behavior efficiently.
The Caputo-Fabrizio derivative is renowned for addressing nonlinear differential equations
efficiently, reducing the computational complexity and the time taken to find the solution
at times. Its solid theoretical foundation and applicability are the keys to popularizing
positive applications of fractional calculus and therefore play a more important role in an
increasingly broad area of science and engineering [22–25] .

Since Rossler has been trailblazing on hyperchaotic systems, there has been significant
advancement to identify multiple applications. Significance includes 2D systems in image
encryption, pseudo-random number generation, and secure communications [26], and [27]
operated in 3D in secure transmission and multi-image encryption, while [28, 29] took
it to the 4D with dynamic analysis and synchronization added. Hyperchaotic systems
have been widely applied in information processing, electronics, neuroscience, and secure
communications, including the encryption of images, audio files, and videos, and also in
the generation of random numbers [30–33].

Fractional-order systems give a better description of complex dynamics, capturing
memory effects and long-term dependencies beyond the reach of traditional models. Frac-



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 3 of 19

tional system solutions improve the simulation accuracy in most fields, and comparative
studies [34, 35] establish the advantages and limitations of different approaches to en-
hance more efficient modeling techniques. Further, chaos analysis of fractional systems
is crucial for safe communication, encryption, and control systems. Modern research on
hyperchaotic systems, which are sensitive and have complicated behavior, is important
for progressing secure data transmission, cryptography, and signal processing. It is also
necessary to study these systems to advance science and technology [36–39].

Fu et al. [40] introduced the following 4-dimensional hyperchaotic system:

dx1
dt

= a(x2 − x1) + x4,

dx2
dt

= cx2 − 10x1x3,

dx3
dt

= −bx3 + 10x1x2,

dx4
dt

= dx2 + x21,

(1)

where a, b, and c are constants, d > 0 is a variable parameter, and x1, x2, x3, and x4 are
driving variables.

Fractional hyperchaotic systems utilizing the Caputo derivative are a generalization of
conventional hyperchaotic systems, described by incorporating fractional-order dynamics
described by: 

CDν
0,tx1 = a(x2 − x1) + x4,

CDν
0,tx2 = cx2 − 10x1x3,

CDν
0,tx3 = −bx3 + 10x1x2,

CDν
0,tx4 = dx2 + x21,

(2)

where CDν
0,t is the Caputo fractional derivative of order ν, while a, b, and c are constants,

d > 0 is a variable parameter, and (x1), (x2), (x3), and (x4) are driving variables.
Despite extensive studies on chaotic dynamics in systems characterized by fractional

orders, there is a lack of appropriate approaches that provide both numerical accuracy and
analytical stability. Many existing approaches are either just numerical simulations lacking
theoretical insight or exclusively analytical approximations that may not sufficiently rep-
resent chaotic behavior. This work fills this gap by integrating the merits of both RPSM
and CFD to provide a well-balanced technique for the analysis of fractional-order chaotic
systems. Moreover, previous research has concentrated much on low-dimensional systems,
while the present study investigates higher-dimensional studies of the 4D fractional-order
system to provide a more precise and realistic approximation of chaotic dynamics. The
results demonstrate that CFD- and RPSM-based methods have a significant potential for
enhancing modeling and solving complicated dynamical systems, indicating their applica-
bility in various scientific and engineering fields [41, 42].



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 4 of 19

This article presents a novel approach, which combines the Caputo fractional deriva-
tive (CFD) technique with the residual power series method (RPSM), to investigate the
complex dynamics within a fractional-order 4D system. Fractional-order chaotic systems
have been studied individually through numerical or analytical methods in the literature
before; however, combined usage of CFD and RPSM in this article presents a more pro-
found and precise insight into chaotic dynamics. The CFD approach guarantees numerical
simulations of high accuracy exactly with numerical certainty, elegantly approximating
complicated chaotic dynamics, while the RPSM guarantees computationally efficient and
stable means of determining analytical solutions. The combination enhances the reliabil-
ity and accuracy of fractional-order chaotic systems, marking a significant advancement
in the field. The paper demonstrates how the two approaches complement each other in
their ways to close the gap between numerical and analytical methods to fractional-order
chaos.

2. Mathematical Preliminaries

Definition 1[43]: The Riemann–Liouville (R-L) fractional integral of a function Q :
R+ → R of order ν > 0 is characterized as follows:

0J
ν
t (Q(t)) =

1

Γ(ν)

∫ t

0
(t− µ)ν−1Q(µ)dµ, t > 0,

where, Γ(ν) represents the Euler gamma function.
Definition 2[43]: The R–L fractional differential operator applied to Q(t) of order

ν > 0 is written as:

RLDν
0,t(Q(t)) =

1

Γ(m− ν)

dm

dtm

∫ t

0
(t− µ)m−ν−1Q(µ)dµ, t > 0,

where, m − 1 < ν < m, andm ∈ N. In this definition, the associated kernel may exhibit
a singularity at one endpoint of the integration interval.

Definition 3[43]: The fractional derivative of a function Q(t) is defined according to
the Caputo approach of a specific order ν > 0 is defined as:

CDν
0,t(Q(t)) =

1

Γ(m− ν)

∫ t

0
(t− µ)m−ν−1Q(m)(µ)dµ, t > 0,

where, m− 1 < ν ≤ m, andm ∈ N.
This Caputo definition is tailored to tackle the difficulties associated with substituting

fractional-order initial conditions using the R-L fractional differential operator in practical
applications. This is accomplished by implementing standard integer-order initial condi-
tions, thus removing the singular kernel found in the R-L definition.
Definition 4: [44] A power series representation of the type

∞∑
m=0

dm(t− t0)
mν = d0 + d1(t− t0)

ν + d2(t− t0)
2ν + ..., (3)



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 5 of 19

where, m − 1 < ν ≤ m and t ≥ t0 is called fractional PS about t0, where t denotes a
variable and dm represents constants known as the coefficients of the series.

Theorem 1: [44] Assume that the function Q possesses a fractional PS representation
at t = t0 characterized by the following form:

Q(t) =
∞∑

m=0

dm(t− t0)
mν , (4)

where, m− 1 < ν ≤ m and t0 ≤ t < t0 +R. If the derivatives DmνQ(t) are continuous in
the interval (t0, t0 + R), where m = 0, 1, 2, ..., then the coefficients dm of Eq. (4) can be
determined using the following formula:

dm =
DmνQ(t0)

Γ(mν + 1)
, (5)

where, Dmν = DνDν · · ·Dν (m-times) and R denotes the radius of convergence.
Definition 5: [44] For m− 1 < ν ≤ m, a power series representation of the type

∞∑
m=0

Qm(x)(t− t0)
mν = Q0(x) +Q1(x)(t− t0)

ν +Q2(x)(t− t0)
2ν + ..., (6)

this is referred to as a multiple fractional PS centered at t = t0, where t represents a
variable and Qm(x) denotes the functions of x which are known as the coefficients of the
series.

3. Numerical Scheme for the Fractional-Order 4D System with Caputo
Derivative

In this part, we present an approximate solution for the model under the fractional
Caputo–Fabrizio operator employing a highly effective numerical approach [45] . The
governing equations are given as:

∗Dν
0,tQ(t) = F (Q(t)), t ∈ [0, a], Q(0) = Q0, (7)

where, Q = (x1, y, z) ∈ R3
+, and the function F (Q) meets the Lipschitz condition:

∥F (Q1(t))− F (Q2(t))∥ ≤ K∥Q1(t)−Q2(t)∥, (8)

with K > 0 as a Lipschitz constant. From Equation (7), we can express the solution
as:

Q(t) = Q0 + 0J
ν
t F (Q(t)), t ∈ [0, a], (9)

where, 0J
ν
t is the fractional integral operator derived from the Caputo operators.

Consider the interval [0, a] divided into n equal parts with step length ∆t = a
n (e.g.,

∆t = 0.05 in simulations). Let x1r be the approximation of x1(t) in t = tr, where r =



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 6 of 19

0, 1, . . . , n. Using finite differences, the numerical scheme for Equation (7) is formulated
as:

cx1r+1
= x10 +

(∆t)ω

Γ(ω + 1)

r∑
k=0

[(1 + r − k)ω − (r − k)ω]F (x1k) +O(∆t2), (10)

Where, O(∆t2) is the order of approximation accuracy.

4. Impact of Fractional Order on the Dynamics of Chaotic Systems via
CFD

In this part, we analyze how the parameter ν influences the behavior of the system, as
illustrated in the figures 1, 2, and 3. Using the Caputo fractional derivative, plots show
how the system changes as the fractional order parameter α takes on different values,
showing that it does have an effect on how the system behaves. From Figure 1, for ν = 1,
the system provides classical chaotic dynamics for integer-order derivatives with densely
structured trajectories in the phase space and the states (x1, x2, x3) exhibiting strong
oscillation activity but no fading out in the time series. Lowering ν to 0.99 in Figure 2;
adds a fractional-order component, hence lowering the system’s chaos with more organized
phase space trajectories. The time series plots have little damping of the oscillations,
showing early energy dissipation. An even smaller attractor is formed by further lowering
alpha to 0.95 in Figure 3; this shows a significant reduction in chaotic behavior. Where the
oscillations of x1, x2, x3 decrease with time significantly, showing a trend towards stability,
the damping is evident even in the time-series plots. Overall, the plots show that lowering
ν shifts the system from a state characterized by classical chaos to a less chaotic and
stable one. This shows that fractional derivatives can be used to stabilize a system and its
behavior. This analysis is important for explaining flow behavior from idealized stability
to actual turbulence in the real world and for controlling the management of systems in
engineering applications.

Figure 1: System dynamics for ν = 1 employing the Caputo fractional derivative.



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 7 of 19

Figure 2: System dynamics for ν = 0.99 employing the Caputo fractional derivative.

Figure 3: System dynamics for ν = 0.95 using the Caputo fractional derivative.

The numerical solutions displayed in Tables 1, 2, and 3 demonstrate the behavior of
the system 2 for various fractional orders (ν = 1, 0.99, 0.95) in t = 2. When ν = 1, the
system functions as an integer-order differential equation, and the CFD method produces
solutions that converge to the 4th-order Runge-Kutta (RK4) reference values with decreas-
ing step size h. As ν decreases to 0.99 and 0.95, the influence of fractional order effects
results in memory dependence, causing a minor reduction in the values of x1, x2, x3, x4.
The CFD method exhibits consistent convergence in all cases, with solutions aligning
closely with the Adams-Bashforth-Moulton method (ABM) reference values for fractional
orders. The findings show that a decrease in h is directly proportional to increased ac-
curacy, which supports the consistency of CFD in solving fractional order systems. The
nature of solutions with decreasing ν represents the effect of fractional dynamics, which
retards the evolution of the system relative to the classical case. The study demonstrates
the efficiency of CFD in the management of systems of integer and fractional order with
high precision.



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 8 of 19

Table 1: Solutions of system 2 by using CFD for (x10, x20, x30, x40) = (1, 1, 1, 1) where
ν = 1, and t = 2

h x1 x2 x3 x4

1/320 0.973059863422542 0.800578111799992 1.591562295686195 2.536897475278398
1/640 0.970285557191089 0.796036503750927 1.593731307460380 2.546833743002026
1/1280 0.969355303143915 0.794508406142886 1.594465238932689 2.550156263728269
1/2560 0.968888986811479 0.793741513533585 1.594834214830244 2.551819454601515
1/5120 0.968608795575532 0.793280445317936 1.595056238342953 2.552817983348334
1/10240 0.967953794109018 0.792201844098075 1.595576140721881 2.555149666256854
RK4 0.967484855614437 0.791428991073163 1.595949065222047 2.556816674404389

Table 2: Solutions of system 2 by using CFD for (x10, x20, x30, x40) = (1, 1, 1, 1) where
ν = 0.99, and t = 2

h x1 x2 x3 x4

1/320 0.958724904468643 0.784005868539537 1.586882323107994 2.552704617046767
1/640 0.957783855967246 0.782526656938586 1.587511271988065 2.555919465717640
1/1280 0.957312261656128 0.781784430216698 1.587827659031779 2.557528691550405
1/2560 0.0.957028942786839 0.781338234845598 1.588018096550532 2.558494798445208
1/5120 0.956839907974789 0.781040408427304 1.588145305517077 2.559139106587635
1/10240 0.955883492764598 0.779514856699298 1.588816353323760 2.562396671913593
ABM 0.955344465148585 0.779339094329801 1.589131949997077 2.567697787348436

Table 3: Solutions of system 2 by using CFD for (x10, x20, x30, x40) = (1, 1, 1, 1) where
ν = 0.95, and t = 2

h x1 x2 x3 x4

1/320 0.915581297407851 0.741241155418083 1.557139066272577 2.575478404944152
1/640 0.915134211844442 0.740621399693080 1.557299568921548 2.576919332084313
1/1280 0.914865705487542 0.740248850488766 1.557396449608853 2.577784401113457
1/2560 0.914686591456969 0.740000191391083 1.557461275589489 2.578361325295418
1/5120 0.913702159352431 0.738585624017066 1.557873715167465 2.581343575041200
1/10240 0.914387867668505 0.739585234231689 1.557569739937361 2.579323238417623
ABM 0.914496444818986 0.739349393578095 1.557768565532089 2.579413747913888

5. Residual power series approach for hyperchaotic system (RPSM)

The fractional hyperchaotic system is recognized as a significant universal nonlinear
model that appears in various physical systems. This article aims to derive the solution
for this equation employing the advantageous RPS method.



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 9 of 19



CDν
0,tx1 = a(x2 − x1) + x4,

CDν
0,tx2 = cx2 − 10x1x3,

CDν
0,tx3 = −bx3 + 10x1x2,

CDν
0,tx4 = dx2 + x21,

(11)

subject to the following constraints initial conditions:

x1(0) = x10, x2(0) = x20, x3(0) = x30, x4(0) = x40. (12)

To numerically address the fractional order 4D hyperchaotic system using the Caputo
derivative, we adhere to the following procedure outlined below.
The RPS method involves representing the solutions of the Eq. (11), given the i.c outlined
in Eq. (12), as a series of various fractional PS expansions centered around the starting
point t = 0. Assume a solution expressed as a power series;

x1(t) =
∞∑
k=0

a1k
tkν

Γ(1 + kν)
,

x2(t) =

∞∑
k=0

a2k
tkν

Γ(1 + kν)
,

x3(t) =

∞∑
k=0

a3k
tkν

Γ(1 + kν)
,

x4(t) =

∞∑
k=0

a4k
tkν

Γ(kν + 1)
.

(13)

The RPS offers analytical approximate solutions expressed as an infinite multiple fractional
PS. To derive the numerical representation of the values depicted in this series, it is
necessary to truncate the resulting series and follow a practical procedure to achieve this
objective. So; let x1n(t), x2n(t), x3n(t), x4n(t) to denot the n-th truncated series solution
of x1(t), x2(t), x3(t), x4(t), in that order. That is :

x1n(t) = x10 +

∞∑
k=0

a1k
tkν

Γ(1 + kν)
,

x2n(t) = x20 +
∞∑
k=0

a2k
tkν

Γ(1 + kν)
,

x3n(t) = x30 +
∞∑
k=0

a3k
tkν

Γ(1 + kν)
,

x4n(t) = x40 +
∞∑
k=0

a4k
tkν

Γ(1 + kν)
.

(14)



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 10 of 19

Define the residual functions for the system :
Res x1(t) =

C Dν
0,tx1 − a(x2 − x1)− x4,

Res x2(t) =
CDν

0,tx2 − cx2 + 10x1x3,

Res x3(t) =
CDν

0,tx3 + bx3 − 10x1x2,

Res x4(t) =
CDν

0,tx4 − dx2 − x21.

(15)

Hence, the n-th residual functions of x1, x2, x3, x4 are :
Res x1n(t) =

C Dν
0,tx1n − a(x2n − x1n)− x4n,

Res x2n(t) =
CDν

0,tx2n − cx2n + 10x1nx3n,

Res x3n(t) =
CDν

0,tx3n + bx3n − 10x1nx2n,

Res x4n(t) =
CDν

0,tx4n − dx2n − x21n.

(16)

Obviously; Res x1(t) = Res x2(t) = Res x3(t) = Res x4(t) = 0, ∀t ≥ 0. So,
limn→∞Res x1n(t) = Res x1(t); limn→∞Res x2n(t) = Res x2(t); limn→∞Res x3n(t) =
Res x3(t); limn→∞Res x4n(t) = Res x4(t). As, the Caputo derivative of any constant is
zero, then : 

CD
(k−1)ν
0,t Res x1(0) =

C D
(k−1)ν
0,t Res x1k(0),

CD
(k−1)ν
0,t Res x2(0) =

C D
(k−1)ν
0,t Res x2k(0),

CD
(k−1)ν
0,t Res x3(0) =

C D
(k−1)ν
0,t Res x3k(0),

CD
(k−1)ν
0,t Res x4(0) =

C D
(k−1)ν
0,t Res x4k(0),

(17)

for k = 1, ..., n.

Now, to obtain the coefficients a1k, a2k, a3k, a4k, k = 1, 2, 3, ..., n, Substituting the n-
truncated series of x1(t), x2(t), x3(t), x4(t) into (13), and then apply the Caputo operators
CD

(n−1)ν
0,t on Res x1n(t), Res x2n(t), Res x3n(t), Res x4n(t) resp., we get :

CD
(n−1)ν
0,t Res x1n(0) = 0,

CD
(n−1)ν
0,t Res x2n(0) = 0,

CD
(n−1)ν
0,t Res x3n(0) = 0,

CD
(n−1)ν
0,t Res x4n(0) = 0,

(18)

for n = 1, 2, 3, ....

Let us find the first few coefficients :

• For n = 1 :



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 11 of 19

x1(t) = x10 + a11
tν

Γ(ν + 1)
= 1 +

a11
Γ(ν + 1)

tν ,

x2(t) = x20 + a21
tν

Γ(ν + 1)
= 1 +

a21
Γ(ν + 1)

tν ,

x3(t) = x30 + a31
tν

Γ(ν + 1)
= 1 +

a31
Γ(ν + 1)

tν ,

x4(t) = x40 + a41
tν

Γ(ν + 1)
= 1 +

a41
Γ(ν + 1)

tν .

(19)

According to (16), the first residual functions of x1(t), x2(t), x3(t), x4(t) are :
Res x11(t) =

C Dν
0,tx11 − a(x21 − x11)− x41,

= a11 − 1− a41t
ν

Γ(ν + 1)
− a

a21t
ν

Γ(ν + 1)
+ a

a11t
ν

Γ(ν + 1)
,

(20)

and
Res x21(t) =

C Dν
0,tx21 − cx2n + 10x1nx3n

= a21 − 2 + d
a11t

ν

Γ(ν + 1)
− c

a21t
ν

Γ(ν + 1)
+ d

a31t
ν

Γ(ν + 1)
+ d

a11a31t
2ν

(Γ(ν + 1))2
,

(21)

also
Res x31(t) =

C Dν
0,tx31 + bx31 − 10x11x21,

= a31 − 7 + b
a31t

ν

Γ(ν + 1)
− d

a21t
ν

Γ(ν + 1)
− d

a11t
ν

Γ(ν + 1)
− d

a11a21t
2ν

(Γ(ν + 1))2
,

(22)

and 
Res x41(t) =

C Dν
0,tx41 − dx21 − x211,

= a41 − 11− d
a21t

ν

Γ(ν + 1)
− 2

a11t
ν

Γ(ν + 1)
− a211t

2ν

(Γ(ν + 1))2
.

(23)

Using (18), we get a11, a21, a31, a41, so that the first RPS solution of Eqs. (11-12) can
be articulated as : 

x11(t) = 1 +
tν

Γ(ν + 1)
,

x21(t) = 1 +
2tν

Γ(ν + 1)
,

x31(t) = 1 +
7tν

Γ(ν + 1)
,

x41(t) = 1 +
11tν

Γ(ν + 1)
.

(24)

• For n = 2 :



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 12 of 19

x12(t) = 1 +
tν

Γ(ν + 1)
+

a12t
2ν

Γ(2ν + 1)
,

x22(t) = 1 +
2tν

Γ(ν + 1)
+

a22t
2ν

Γ(2ν + 1)
,

x32(t) = 1 +
7tν

Γ(ν + 1)
+

a32t
2ν

Γ(2ν + 1)
,

x42(t) = 1 +
11tν

Γ(ν + 1)
+

a42t
2ν

Γ(2ν + 1)
,

(25)

and the second residual functions are :
Res x12(t) =

C Dν
0,tx12 − a(x22 − x12)− x42,

= (a22−46)
tν

Γ(ν + 1)
− (a42 + a(a22 − a12)))

t2ν

Γ(2ν + 1)
,

(26)

and

Res x22(t) =
C Dν

0,tx22 − cx22 + dx12x32

= (a22 + 56)
tν

Γ(ν + 1)
+ (d(a32 + a12 − 12a22))

t2ν

Γ(2ν + 1)
+ 70

t2ν

(Γ(ν + 1))2

+ (da32 + 70a12)
t3ν

Γ(ν + 1)Γ(2ν + 1)
+

da12a32t
4ν

(Γ(2ν + 1))2
,

(27)
also

Res x32(t) =
C Dν

0,tx32 + bx32 − 10x12x22,

= (a32 − 9)
tν

Γ(ν + 1)
+ (ba32 − d(a22 + a12))

t2ν

Γ(2ν + 1)
− 20

t2ν

(Γ(ν + 1))2

− (da22 + 20a12)
t3ν

Γ(ν + 1)Γ(2ν + 1)
− da12a22

t4ν

(Γ(2ν + 1))2
,

(28)
and

Res x42(t) =
C Dν

0,tx42 − dx22 − x212,

= (a42 − 22)
tν

Γ(ν + 1)
− (d+ 2a12

t2ν

Γ(2ν + 1)
− t2ν

(Γ(ν + 1))2
− 2a12

t3ν

Γ(ν + 1)Γ(2ν + 1)

− a212t
4ν

(Γ(2ν + 1))2
.

(29)

Using (18), we get a12, a22, a32, a42, therefore, the second RPS solution of Eqs. (11-
12) are given by the following form :



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 13 of 19



x12(t) = 1 +
tν

Γ(ν + 1)
+ 46

t2ν

Γ(2ν + 1)
,

x22(t) = 1 +
2tν

Γ(ν + 1)
− 56

t2ν

Γ(2ν + 1)
,

x32(t) = 1 +
7tν

Γ(ν + 1)
+ 9

t2ν

Γ(2ν + 1)
,

x42(t) = 1 +
11tν

Γ(ν + 1)
+ 22

t2ν

Γ(2ν + 1)
.

(30)

• In the same manner, implementing the identical procedures for n = 3, leads to :

x13(t) = 1 +
tν

Γ(ν + 1)
+ 46

t2ν

Γ(2ν + 1)
+ a13

t3ν

Γ(3ν + 1)
,

x23(t) = 1 +
2tν

Γ(ν + 1)
− 56

t2ν

Γ(2ν + 1)
+ a23

t3ν

Γ(3ν + 1)
,

x33(t) = 1 +
7tν

Γ(ν + 1)
+ 9

t2ν

Γ(2ν + 1)
+ a33

t3ν

Γ(3ν + 1)
,

x43(t) = 1 +
11tν

Γ(ν + 1)
+ 22

t2ν

Γ(2ν + 1)
+ a43

t3ν

Γ(3ν + 1)
.

(31)

By the same process, we get :

a13 = 3548, a23 = −1222− 70Γ(2ν+1)
(Γ(ν+1))2

, a33 = −127+ 20Γ(2ν+1)
(Γ(ν+1))2

, a43 = −514+ Γ(2ν+1)
(Γ(ν+1))2

.

Hence,

x13(t) = 1 +
tν

Γ(ν + 1)
+ 46

t2ν

Γ(2ν + 1)
+ 3548

t3ν

Γ(3ν + 1)
,

x23(t) = 1 +
2tν

Γ(ν + 1)
− 56

t2ν

Γ(2ν + 1)
− (1222 +

70Γ(2ν + 1)

(Γ(ν + 1))2
)

t3ν

Γ(3ν + 1)
,

x33(t) = 1 +
7tν

Γ(ν + 1)
+ 9

t2ν

Γ(2ν + 1)
+ (−127 +

20Γ(2ν + 1)

(Γ(ν + 1))2
)

t3ν

Γ(3ν + 1)
,

x43(t) = 1 +
11tν

Γ(ν + 1)
+ 22

t2ν

Γ(2ν + 1)
+ (−514 +

Γ(2ν + 1)

(Γ(ν + 1))2
)

t3ν

Γ(3ν + 1)
.

(32)

This process can be reiterated until the coefficients associated with the multiple frac-
tional PS solutions of equations (11-12) are acquired in the desired order.

Here, the parameter values are taken as a = 35, b = 3, c = 12, and d = 10.



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 14 of 19

6. Application of the CFD and RPSM to Chaotic Dynamical System

In this section, we present numerical solutions acquired with the Residual Power Series
Method (RPSM) to solve the system 2. The solutions are computed for different values
of fractional order (0.95, 0.99, and 1) with homogeneous initial conditions. The objec-
tive is to investigate the effect of fractional order on the dynamics of the system and to
analyze the accuracy and convergence of the RPSM method. Tables 4, 5, and 6 exhibit
numerical data for the first 4 terms of the solutions, showing solutions of the system 2
employing the Residual Power Series Method (RPSM) for various fractional order values
of ν (0.95, 0.99, and 1). The data show a decrease with decreasing t values and eventu-
ally converge, hence verifying stability and the correctness of the approach. Different ν
values compared reveal little variance in the solutions, which result from the fractional-
order parameter on system behavior. The proximity of values with diminishing step sizes
further attests to the effectiveness of RPSM. These results confirm the method’s efficacy
in fractional-order differential systems, hence proving its relevance in mathematical mod-
eling.

Table 4: Solutions of system (2) using RPSM for (x10, x20, x30, x40) = (1, 1, 1, 1) where
ν = 1.

t x1 x2 x3 x4

0.5 0.501189960570766 1.010000344410005 3.540065400660771 5.564300000813071
0.6 0.600000699400806 1.210000155410008 4.200061400760091 6.663400200643092
0.7 0.700460199106844 1.441241155418083 4.966107145270071 7.764400204643071
0.8 0.821342118443442 1.630061378692061 5.665176578821620 8.885010322064221
0.9 0.933755705487631 1.850027650460755 6.369196439800855 9.976584402223558
1.0 1.009700731466999 2.000800191391083 7.000761275589489 11.00897725296401

Table 5: Solutions of system (2) using RPSM for (x10, x20, x30, x40) = (1, 1, 1, 1) where
ν = 0.99.

t x1 x2 x3 x4

0.5 0.505595555300738 1.011193002380041 3.539167400662136 5.561546530022142
0.6 0.605608366300064 1.211224063300052 4.239265004152555 6.661698300312111
0.7 0.705455377406621 1.410910065307072 4.938186043162656 7.761698403922141
0.8 0.805157220074431 1.809471378646070 6.333169556721431 9.952118122074211
0.9 0.904737604366541 1.809476670465756 6.333165237606752 9.952162304313457
1.0 1.004275561446867 2.008410162371001 7.029432265578476 11.04625221527541



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 15 of 19

Table 6: Solutions of system (2) using RPSM for (x10, x20, x30, x40) = (1, 1, 1, 1) where
ν = 0.95.

t x1 x2 x3 x4

0.5 0.528261457200942 1.056520060720154 3.697832000373467 5.810878940022073
0.6 0.628160365000042 1.256321072320061 4.397125003153668 6.909767440722251
0.7 0.727227431023701 1.454451329793310 5.090599567821537 7.999519332062421
0.8 0.825854310237016 1.651171388672142 5.779199568921548 9.081457211094202
0.9 0.923218096487542 1.454457750267755 5.090595337506652 7.999574401223463
1.0 1.020536481457856 2.041060182371174 7.143732264476477 11.22595132227552

7. Conclusion

To investigate the intricate dynamics within the four-dimensional fractional order sys-
tem, this study shows the efficiency of the residual power series method (RPSM) and the
Caputo fractional derivative (CFD). The RPSM is a reliable and simple way to get analyti-
cal answers; the CFD approach is seen to provide extremely accurate numerical simulations
and effectively capture chaotic behavior, and by means of these two approaches, one gains
additional insight into the complex dynamics of the system. In addition to providing use-
ful tools for modeling and solving fractional-order systems, the findings reveal that both
techniques are successful in discovering chaotic attractors and closely approximating the
behavior of these systems. Due to the accuracy and consistency of these methodologies,
they offer tremendous potential for future usage in engineering and scientific research,
particularly in areas of study where fractional-order models play a key role. In the future,
we intend to solve some new fractional models, such as in [46–48] and make comparisons
with other numerical methods [49–51].

Acknowledgements

The authors thank the deanship in Zarqa University. This research is fully funded by
Zarqa University-Jordan.

Conflicts of Interest: The authors declare that they have no conflict of interest.

References

[1] Yong Zhang et al. A review of applications of fractional calculus in earth system
dynamics. Chaos, Solitons & Fractals, 102:29–46, 2017.

[2] Yuriy A. Rossikhin and Marina V. Shitikova. Application of fractional calculus for
dynamic problems of solid mechanics: Novel trends and recent results. Applied Me-
chanics Reviews, page 010801, 2010.

[3] Ivo Petras. Fractional calculus and its applications. Mathematical Modeling with
Multidisciplinary Applications, pages 355–396, 2012.



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 16 of 19

[4] Vasily E. Tarasov. On history of mathematical economics: Application of fractional
calculus. Mathematics, 7(6):509, 2019.

[5] Faeza Hasan et al. A new perspective on the stochastic fractional order materialized
by the exact solutions of allen-cahn equation. International Journal of Mathematical,
Engineering and Management Sciences, 8(5):912, 2023.

[6] Sana Abdulkream, Alharbi et al. Modeling and analysis of visceral leishmaniasis dy-
namics using fractional-order operators: A comparative study. Mathematical Methods
in the Applied Sciences, 47(12):9918–9937, 2024.

[7] Rania Saadeh et al. Mathematical modeling and stability analysis of the novel frac-
tional model in the caputo derivative operator: A case study. Heliyon, 10(5), 2024.

[8] Faeza Lafta Hasan et al. Exploring analytical results for (2+1) dimensional breaking
soliton equation and stochastic fractional broer-kaup system. AIMS Mathematics,
9(5):11622–11643, 2024.

[9] Dalal Khalid Almutairi et al. A numerical confirmation of a fractional seitr for in-
fluenza model efficiency. Applied Mathematics, 17(5):741–749, 2023.

[10] Mohamed A. Abdoon and Abdulrahman B. M. Alzahrani. Comparative analysis
of influenza modeling using novel fractional operators with real data. Symmetry,
16(9):1126, 2024.

[11] E. Nahaa Alsubaie et al. Improving influenza epidemiological models under caputo
fractional-order calculus. Symmetry, 16(7):929, 2024.

[12] Dalal Khalid Almutairi et al. Modeling and analysis of a fractional visceral leishman-
iosis with caputo and caputo-fabrizio derivatives. Journal of the Nigerian Society of
Physical Sciences, pages 1453–1453, 2023.

[13] Mohamed Elbadri et al. A numerical solution and comparative study of the symmetric
rossler attractor with the generalized caputo fractional derivative via two different
methods. Mathematics, 11(13):2997, 2023.

[14] Sérgio Adriani David, Juan Lopez Linares, and Eliria Maria de Jesus Agnolon Pallone.
Fractional order calculus: Historical apologia, basic concepts and some applications.
Revista Brasileira de Ensino de F́ısica, 33:4302–4302, 2011.

[15] Abdon Atangana and J. F. Gómez-Aguilar. Numerical approximation of riemann-
liouville definition of fractional derivative: From riemann-liouville to atangana-
baleanu. Numerical Methods for Partial Differential Equations, 34(5):1502–1523,
2018.

[16] Obaid Jefain Julaighim Algahtani. Comparing the atangana-baleanu and caputo-
fabrizio derivative with fractional order: Allen cahn model. Chaos, Solitons & Frac-
tals, 89:552–559, 2016.

[17] Nadeem Ahmad Sheikh et al. A comparative study of atangana-baleanu and caputo-
fabrizio fractional derivatives to the convective flow of a generalized casson fluid. The
European Physical Journal Plus, 132:1–14, 2017.

[18] Adil Khurshaid and Hajra Khurshaid. Comparative analysis and definitions of frac-
tional derivatives. Journal ISSN, 2766:2276, 2023.

[19] A. Chavada and N. Pathak. Transmission dynamics of breast cancer through caputo
fabrizio fractional derivative operator. Mathematical Modelling and Control, 4:119–



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 17 of 19

132, 2024.
[20] Abdon Atangana and Kolade M. Owolabi. New numerical approach for fractional

differential equations. Mathematical Modelling of Natural Phenomena, 13(1):3, 2018.
[21] Ilhem Kadri, Mohammed Al Horani, and Roshdi Khalil. Solution of fractional laplace

type equation in conformable sense using fractional fourier series with separation of
variables technique. Results in Nonlinear Analysis, 6(2):53–59, 2023.

[22] S. Abdulkream Alharbi, M. A. Abdoon, R. Saadeh, R. Allogmany, M. Berir, and
F. EL Guma. Modeling and analysis of visceral leishmaniasis dynamics using
fractional-order operators: A comparative study. Mathematical Methods in the Ap-
plied Sciences, 47(12):9918–9937, 2024.

[23] A. B. M. Alzahrani, R. Saadeh, M. A. Abdoon, M. Elbadri, M. Berir, and A. Qazza.
Effective methods for numerical analysis of the simplest chaotic circuit model with
atangana-baleanu caputo fractional derivative. Journal of Engineering Mathematics,
144(1):9, 2024.

[24] A. Qazza and R. Saadeh. On the analytical solution of fractional sir epidemic model.
Applied Computational Intelligence and Soft Computing, 2023:16, 2023.

[25] E. Salah, A. Qazza, R. Saadeh, and A. El-Ajou. A hybrid analytical technique for
solving multi-dimensional time-fractional navier-stokes system. AIMS Mathematics,
8(1):1713–1736, 2023.

[26] Uğur Erkan, Abdurrahim Toktas, and Qiang Lai. 2d hyperchaotic system based on
schaffer function for image encryption. Expert Systems with Applications, 213:119076,
2023.

[27] A. Sahasrabuddhe and D. S. Laiphrakpam. Multiple images encryption based on 3d
scrambling and hyper-chaotic system. Information Sciences, 550:252–267, 2021.

[28] H. Wang and G. Dong. New dynamics coined in a 4-d quadratic autonomous hyper-
chaotic system. Applied Mathematics and Computation, 346:272–286, 2019.

[29] T. Nestor, A. Belazi, B. Abd-El-Atty, M. N. Aslam, C. Volos, N. J. De Dieu, et al. A
new 4d hyperchaotic system with dynamics analysis, synchronization, and application
to image encryption. Symmetry, 14:424, 2022.

[30] ShiMing Fu, XueFeng Cheng, and Juan Liu. Dynamics, circuit design, feedback con-
trol of a new hyperchaotic system and its application in audio encryption. Scientific
Reports, 13(1):19385, 2023.

[31] Xiaodong Li et al. Video encryption based on hyperchaotic system. Multimedia Tools
and Applications, 79:23995–24011, 2020.

[32] A. Hadj Brahim et al. A novel pseudo-random number generator: Combining hyper-
chaotic system and des algorithm for secure applications. The Journal of Supercom-
puting, 81(1):94, 2025.

[33] Yang Liu and Xiaojun Tong. Hyperchaotic system-based pseudorandom number gen-
erator. IET Information Security, 10(6):433–441, 2016.

[34] Ilhem Kadri, Mohammed Al Horani, and Roshdi R. Khalil. Solution of non-linear
fractional burger’s type equations using the laplace transform decomposition method.
Results in Nonlinear Analysis, 5(2):131–150, 2022.

[35] H. Gündoğdu and H. Joshi. Numerical analysis of time-fractional cancer models with



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 18 of 19

different types of net killing rate. Mathematics, 13(3):536, 2025.
[36] Muhammad Sarfraz et al. Dynamic analysis of a 10-dimensional fractional-order

hyperchaotic system using advanced hyperchaotic metrics. Fractal and Fractional,
9(2):76, 2025.

[37] Thwiba Abdulrhman. Stability analysis of fractional chaotic and fractional-order
hyperchain systems using lyapunov functions. European Journal of Pure and Applied
Mathematics, 18(1):5576–5576, 2025.

[38] Ammar Soukkou et al. Design tools to stabilize and to synchronize fractional-order en-
ergy resources system based on fractional-order control approaches: A review. Journal
of the Brazilian Society of Mechanical Sciences and Engineering, 47(4):173, 2025.

[39] J. Jackson and R. Perumal. A robust image encryption technique based on an im-
proved fractional order chaotic map. Nonlinear Dynamics, 113(7):7277–7296, 2025.

[40] S. Fu, X. F. Cheng, and J. Liu. Dynamics, circuit design, feedback control of a
new hyperchaotic system and its application in audio encryption. Scientific Reports,
13:19385, 2023.

[41] H. Gundogdu and O. F. Gozukizil. On the approximate numerical solutions of frac-
tional heat equation with heat source and heat loss. Thermal Science, 26(5 Part
A):3773–3786, 2022.

[42] H. Gündoğdu and Ö. F. Gözükızıl. Applications of the decomposition methods to
some nonlinear partial differential equations. New Trends in Mathematical Sciences,
6(3):57–66, 2018.

[43] M. Caputo and M. Fabrizio. A new definition of fractional derivative without singular
kernel. Progress in Fractional Differentiation and Applications, 1:73–85, 2015.

[44] Omar Abu Arqub et al. A reliable analytical method for solving higher-order initial
value problems. Discrete Dynamics in Nature and Society, 2013(1):673829, 2013.

[45] Sania Qureshi and Abdon Atangana. Mathematical analysis of dengue fever outbreak
by novel fractional operators with field data. Physica A: Statistical Mechanics and
its Applications, 526:121127, 2019.

[46] R. Saadeh, A. Alshawabkeh, R. Khalil, M. A. Abdoon, N. Taha, and D. K. Almutairi.
The mohanad transforms and their applications for solving systems of differential
equations. European Journal of Pure and Applied Mathematics, 17(1):385–409, 2024.

[47] F. E. L. Gumaa, M. A. Abdoon, A. Qazza, R. Saadeh, M. A. Arishi, and A. M. De-
goot. Analyzing the impact of control strategies on visceral leishmaniasis: A mathe-
matical modeling perspective. European Journal of Pure and Applied Mathematics,
17(2):1213–1227, 2024.

[48] R. Allogmany, N. A. Almuallem, R. D. Alsemiry, and M. A. Abdoon. Exploring
chaos in fractional order systems: A study of constant and variable-order dynamics.
Symmetry, 17(4):605, 2025.

[49] M. Abu-Ghuwaleh, R. Saadeh, and A. Qazza. General master theorems of integrals
with applications. Mathematics, 10(19):3547, 2022.

[50] D. E. Elgezouli, H. Eltayeb, and M. A. Abdoon. Novel gpid: Grünwald-letnikov
fractional pid for enhanced adaptive cruise control. Fractal and Fractional, 8(12):751,
2024.



I.Kadri et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6381 19 of 19

[51] M. Ali, S. M. Alzahrani, R. Saadeh, M. A. Abdoon, A. Qazza, N. Al-kuleab, and
F. EL Guma. Modeling covid-19 spread and non-pharmaceutical interventions in
south africa: A stochastic approach. Scientific African, 24:e02155, 2024.