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

Employing the Limit Residual Function Method to
Solve Systems of Fractional Differential Equations

Ahmad El-Ajou1, Aliaa Burqan2,∗

1 Department of Mathematics, Faculty of Science, Al Balqa Applied University, Salt 19117,
Jordan
2 Department of Mathematics, Faculty of Science, Zarqa University, Zarqa 13110, Jordan

Abstract. This study aims to solve systems of fractional differential equations analytically using
a simple new technique, the limit residual function method. This method relies on coupling the
residual function with the limit to produce analytical and approximate solutions within rapidly
converging series forms. This technique could be an alternative to the residual power series method
which is an efficient and quick-to-solve system of fractional differential equations, both linear
and nonlinear, that arise in numerous physical phenomena. To illustrate the methodology and
confirm its effectiveness, the study explores three different applications. The proposed algorithm’s
reliability and accuracy can be easily described by contrasting the numerical results with the exact
solutions.

2020 Mathematics Subject Classifications: 26A33, 32A05, 34A12

Key Words and Phrases: Fractional initial value problems, Fractional power series, Caputo’s
derivative operator, residual function

1. Introduction

In applied sciences such as applied mathematics, mathematical biology, physics, and
engineering, fractional calculus is currently an area of intense study. There are multiple
definitions for the fractional derivative, making it a non-unique formula now. The most
commonly used definition is the Riemann-Liouville (R-L) definition, while Caputo’s def-
inition (1967) of the fractional derivative is also widely used [6, 12, 13, 29, 31, 32, 37].
Differential equations of fractional order have been the subject of several investigations
due to their frequent appearance in different scientific applications. The differential equa-
tions in several forms of fractional derivatives give different types of solutions. As a result,
fractional differential equations cannot be solved using a standard approach. Therefore,
a rapidly developing field of applied mathematics is the interpretation and solution of
fractional differential equations. The Predictor-Corrector approach [16], the variational

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i1.5617

Email addresses: ajou44@bau.edu.jo (A. El-Ajou), aliaaburqan@zu.edu.jo (A. Burqan)

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



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 2 of 19

iteration approach [14, 23], the Homotopy perturbation approach [21], the Homotopy anal-
ysis method [22], the Adomian decomposition method [1, 26, 28], the operational matrix
approach [34], the residual power series method [15, 20, 30], the Laplace residual power
series method [3, 9, 10, 27, 36], the differential transformation technique [11], and the
Taylor series method [25], are newly employed techniques for solving linear and nonlinear
differential equations.

Recently, El-Ajou and Burgan introduced the limit residual function (LRF) method
[18, 19] to find analytical series solutions for linear and nonlinear partial differential equa-
tions. The LRF approach is a semi-analytical technique that constructs an analytical
solution using polynomials based on power series (PS) expansion, with the series coef-
ficients determined by calculating the limit of the residual functions. Indeed, the LRF
offers a power series solution for differential equations that is the same as the Taylor series
solution. LRF is just a new, simple, and efficient technique for determining the coefficients
of the Taylor series.

On the other hand, much literature has developed concerning the fractional system of
differential equations and its applications [1, 2, 4, 5, 26, 35]. In this article, we present
a new application of the LRF technique to yield approximate solutions for the system of
fractional differential equations on the form

Dαu1 (x) = H1 (x, u1 (x) , u2 (x) , . . . , un (x)) ,

Dαu2 (x) = H2 (x, u1 (x) , u2 (x) , . . . , un (x)) ,

...

Dαun (x) = H1 (x, u1 (x) , u2 (x) , . . . , un (x)) , (1)

where Dα, 0 < α ≤ 1 is the derivative of order α in the sense of Caputo, subject to the
initial conditions u1 (0) = u1,0, u2 (0) = u2,0, . . . , un (0) = un,0.

The structure of this paper is as follows: Section 2 revisits several important funda-
mental results of fractional calculus and fractional PS. Section 3 provides an overview of
the suggested method for solving the system of fractional differential equations. In Section
4, three systems of fractional initial value problems are solved to explore the applicability
and simplicity of the LRF technique. Section 5 provides the conclusion.

2. Essential Preliminaries and Notations

This section revisits several important fundamental results about the fractional PS,
essential, to developing analytical solutions for systems of fractional differential equations
using the LRF method. There are many definitions of fractional derivatives in the mathe-
matical literature, such as the Riemann–Liouville fractional derivative, Caputo fractional
derivative, Grünwald–Letnikov fractional derivative [12, 31], the two-scale fractal deriva-
tive [7], He’s fractional derivative [24], conformable fractional derivatives, and Atangana-
Baleanu’s fractional derivative [8]. This article focuses on the Caputo derivative of order
α defined for the function u(x) as in the next definition. In the future, researchers may be



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 3 of 19

able to adapt the LRF method to solve fractional differential equations using other types
of fractional derivatives.

Definition 2.1. [12, 31] The Caputo derivative of order α ∈ (m − 1,m], m∈N of the
function u (x) is defined as follows:

Dαu (x) =

{
1

Γ (α)

∫ x
y0
(x− y)α−1 u(m) (y) dy, m− 1 < α < m , x > y ≥ y0 ≥ 0,

u(m) (x) , α = m.
(2)

Many properties of Caputo derivative can be found in the Refs. [12, 31].

In the following, the definitions and theorems relevant to the classical PS are expanded
to include the fractional case in the Caputo sense.

Definition 2.2. [17] A PS representation of the form

∞∑
j=0

uj (x− x0)
jα = u0 + u1(x− x0)

α + u2 (x− x0)
2α + . . . , (3)

where 0 ≤ m − 1 < α ≤ m, x ≥ x0 is called a fractional PS about x0 where x is variable
and uj ’s are constants called the coefficients of the series.

Theorem 2.1. [33] The fractional PS
∑∞

j=0 uj (x− x0)
jα = 0 for all x in |x− x0| < s if

and only if each coefficient uj equals zero.

To determine the coefficients of the fractional PS solution, the primary principle for
the LRFD approach has been presented and validated by the authors in references [18, 19],
which is depicted in the following theorem.

Theorem 2.2. Suppose that u (x) has the fractional PS expansion u (x) =
∑∞

j=0 uj (x− x0)
jα

and u (x) = 0 for all x in some interval I. Then

lim
x→x0

uk(x)

(x− x0)
(k−1)α

= 0, k = 1, 2, . . . , x ̸= x0, (4)

where

uk(x) =
k∑

j=0

uj (x− x0)
jα. (5)

3. The methodology of the LRF method for solving systems of
fractional differential equations

This section introduces the fundamental idea of the LRF technique for analytically
solving specific systems of linear and non-linear fractional differential equations. The



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 4 of 19

residual function and the limit at zero are the foundational concepts of this approach. To
demonstrate the procedures of the LRF method, we will be studying the following class
of fractional differential equation systems:

Dαui (x) = Hi (x, u1 (x) , u2 (x) , . . . , un (x)) , x ≥ 0, 0 < α ≤ 1, (6)

with the initial conditions

ui (0) = ui,0 , i = 1, 2, . . . , n. (7)

where α is the order of the Caputo fractional differential operator Dα, Hi are analytic
functions, and ui (x) are unknown analytical smooth functions that will be identified.

The LRF method assumes writing the functions ui (x) using the following fractional
PS expansions:

ui (x) =
∞∑
j=0

ui,j xjα, i = 1, 2, . . . , n, x ≥ 0. (8)

Since ui,0 = ui (0) , and by truncating the series (8), we obtain the kth approximate
solution of the system (6)-(7) as follows:

uik (x) = ui,0 +

k∑
j=1

ui,j xjα, i = 1, 2, . . . , n, x ≥ 0, (9)

in which the constants of the expansion series, ui,n, can be obtained by identifying the
residual functions as follows:

Rf (ui (x)) = Dαui (x)−Hi (x, u1 (x) , u2 (x) , . . . , un (x)) , i = 1, 2, . . . , n, x ≥ 0. (10)

So, the kth residual functions can be written as:

Rf (uik (x)) = Dαuik (x)−Hi (x, u1k (x) , u2k (x) , . . . , unk (x)) , i = 1, 2, . . . , n, x ≥ 0.
(11)

To find the kth approximate solution of the system (6)-(7), we need to determine the
value of the coefficients, ui,j , i = 1, 2, . . . , n, and j = 1, 2, . . . , k, in the series (9). The
essential tool of the LRF method, which effectively identifies the unknown coefficients, as
in [25], is

lim
x→0

Rf (uik (x))

x(k−1)α
= 0, i = 1, 2, ..., n. (12)

To find the 1st approximate solution, we must determine the coefficients ui,1 in the
series (9). We can do this by substituting ui1 (x) = ui,0 + ui,1x

α into Rf (ui1 (x)) , i =
1, 2, ..., n, and then solving the equations limx→0Rf (ui1 (x)) = 0 for ui,1, i = 1, 2, . . . , n.
This will lead us to our desire.

Similarly, we can get the actual 2nd approximate solution by substituting ui2 (x) =

ui,0 + ui,1x
α + ui,2x

2α into Rf (ui2 (x)) and solving the equations limx→0
Rf(ui2(x))

xα = 0
for ui,2, i = 1, 2, . . . , n.

Generally, the kth approximate solution of the system (6)-(7) is determined by sub-
stituting uik (x) = ui(k−1) (x) + ui,kx

kα into Rf (uik (x)) and then solving the algebraic
equations (12) for ui,k.



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 5 of 19

4. Illustrative examples

This section tests the performance and implementation of the suggested method on
four systems of linear and non-linear fractional differential equations. The accuracy of the
technique will be evaluated by comparing the results obtained with the exact solutions.

Example 4.1. Consider the following system of linear fractional differential equations:

Dαu1 (x) = u1 (x) + u2 (x) ,

Dαu2 (x) = −u1 (x) + u2 (x) , (13)

where 0 < α ≤ 1, 0 ≤ x with the initial conditions

u1 (0) = 0, u2 (0) = 1. (14)

In the classical case (α = 1), the exact solution may be obtained analytically and
provided as follows:

u1 (x) = exsinx , u2 (x) = excosx . (15)

Considering the LRF technique to create an analytical series solution for System (13)–
(14), we start by assuming that the solution has the following fractional series expansions:

u1 (x) =

∞∑
n=0

u1,n xnα, u2 (x) =

∞∑
n=0

u2,n xnα. (16)

Based on the initial conditions specified in Equation (14), we have u1,0 = 0, u2,0 = 1.
Thus, the kth approximation of u1 (x) and u2 (x) can be expressed as:

u1k (x) =
k∑

n=1

u1,n xnα, u2k (x) = 1 +
k∑

n=1

u2,n xnα. (17)

To identify the additional unknown coefficients of the series given in Equation (16),
we proceed with the second step of the LRF approach by defining the residual functions
of equations in (13) as follows:

Rf (u1 (x)) = Dαu1 (x)− u1 (x)− u2 (x) ,

Rf (u2(x)) = Dαu2 (x) + u1 (x)− u2 (x) , (18)

and the kth residual functions are expressed as follows:

Rf (u1k (x)) = Dαu1k (x)− u1k (x)− u2k (x) ,

Rf (u2k (x)) = Dαu2k (x) + u1k (x)− u2k (x) . (19)



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 6 of 19

By substituting the first approximations u11 (x) = u1,1 xα, u21 (x) = 1 + u2,1x
α into

the first residual functions, Rf (u11 (x)) , Rf (u21 (x)) , we obtain

Rf (u11 (x)) = u1,1Γ (α+ 1)− u1,1x
α − 1− u2,1x

α,

Rf (u21 (x)) = u2,1Γ (α+ 1) + u1,1x
α − 1− u2,1x

α. (20)

Solving the equations limx→0Rf (u11 (x)) = 0, limx→0Rf (u21 (x)) = 0 for u1,1, u2,1,
respectively, we have

u1,1 =
1

Γ (α+ 1)
, u2,1 =

1

Γ (α+ 1)
. (21)

In order to determine the second coefficients u1,2 and u2,2, substitute the second ap-
proximations u12 (x) =

xα

Γ (α+1) +u1,2x
2α, u22 (x) = 1+ xα

Γ (α+1) +u2,2x
2α into Rf (u12 (x)) ,

Rf (u22 (x)) to have

Rf (u12 (x)) = 1 + u1,2
Γ (2α+ 1) xα

Γ (α+ 1)
−
(

xα

Γ (α+ 1)
+ u1,2 x2α

)

−
(
1 +

xα

Γ (α+ 1)
+ u2,2 x2α

)
,

Rf (u22 (x)) = 1 + u2,2
Γ (2α+ 1) xα

Γ (α+ 1)
+

(
xα

Γ (α+ 1)
+ u1,2 x2α

)
−
(
1 +

xα

Γ (α+ 1)
+ u2,2 x2α

)
. (22)

Making simple calculations, the limits limx→0
Rf(u12(x))

xα = 0, limx→0
Rf(u22(x))

xα = 0,
yield

u1,2 =
2

Γ (2α+ 1)
, u2,2 = 0. (23)

The values of the coefficients u1,3 and u2,3 are also obtained by substituting u13 (x) =
xα

Γ (α+1) +
2

Γ (2α+1)x
2α + u1,3 x3α and u23 (x) = 1 + xα

Γ (α+1) + u2,3 x3α into Rf (u13 (x)) ,

Rf (u23 (x)), and then we solve the following equations:

lim
x→0

Rf (u13 (x))

x2α
= u1,3

Γ (3α+ 1)

Γ (2α+ 1)
− 2

Γ (2α+ 1)
= 0,

lim
x→0

Rf (u23 (x))

x2α
= u2,3

Γ (3α+ 1)

Γ (2α+ 1)
+

2

Γ (2α+ 1)
= 0. (24)

Then we have

u1,3 =
2

Γ (3α+ 1)
, u2,3 = − 2

Γ (3α+ 1)
. (25)

Proceeding in the same manner, we get

u1,4 = 0, u2,4 = − 4

Γ (4α+ 1)
. (26)



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 7 of 19

Table 1: The coefficients of the 7th approximation of u1 (x) , u2 (x) for the System (13)-
(14).

k u1,k u2,k

0 0 1

1 1
Γ (α+1)

1
Γ (α+1)

2 2
Γ (2α+1) 0

3 2
Γ (3α+1) − 2

Γ (3α+1)

4 0 − 4
Γ (4α+1)

5 − 4
Γ (5α+1) − 4

Γ (5α+1)

6 − 8
Γ (6α+1) 0

7 − 8
Γ (7α+1)

8
Γ (7α+1)

For k = 5, 6, 7, the algebraic equations limx→0
Rf(u1k(x))

x(k−1)α = 0, limx→0
Rf(u2k(x))

x(k−1)α = 0
can be solved repeatedly to give the coefficients of the 7th approximation. The needed
coefficients are summarized in Table 1.

So, the LRF solution of the System (13)-(14) has the following series expansions:

u1 (x) =
xα

Γ (α+ 1)
+

2x2α

Γ (2α+ 1)
+

2x3α

Γ (3α+ 1)
− 4 x5α

Γ (5α+ 1)
− 8 x6α

Γ (6α+ 1)
− 8 x7α

Γ (7α+ 1)
+. . . ,

u2 (x) = 1+
xα

Γ (α+ 1)
− 2x3α

Γ (3α+ 1)
− 4 x4α

Γ (4α+ 1)
− 4 x5α

Γ (5α+ 1)
+

8 x7α

Γ (7α+ 1)
+ . . . . (27)

When α = 1, the series solution (27) becomes as follows:

u1 (x) = x+ x2 +
x3

3
− x5

30
− x6

90
− x7

630
+ . . . ,

u2 (x) = 1 + x− x3

3
− x4

6
− x5

30
+

x7

630
+ . . . , (28)

which is the expansion of the exact solution given in (15).
In Figure 1, the graphs show the 7th approximate solution of the system (13)-(14) at

different values of α, as well as the exact solution at α = 1. The graph demonstrates
strong agreement between the 7th approximate solution and the exact solution at α = 1.
Furthermore, it illustrates the impact of the derivative order on the solution behavior, re-
vealing that the solution curve decreases as the order of the fractional derivative decreases.



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 8 of 19

(a) (b)

Figure 1: The curves of the 7th approximate and exact (solid curve) solutions of System
(13)-(14) at different values of α.

Example 4.2. Consider the following system of non-linear fractional differential equations

Dαu1 (x) = u1 (x) ,

Dαu2 (x) = 2u21 (x) ,

Dαu3 (x) = 3u1 (x) u2 (x) , (29)

where 0 < α ≤ 1, 0 ≤ x with the initial conditions

u1 (0) = 1, u2 (0) = 1, u3 (0) = 1. (30)

In the classical case (α = 1), the exact solution may be obtained analytically and
provided as follows:

u1 (x) = ex, u2 (x) = e2x, u3 (x) = e3x . (31)

Assume the solution of the System (29)-(30) has the following fractional series expan-
sions:

u1 (x) =

∞∑
n=0

u1,n xnα, u2 (x) =

∞∑
n=0

u2,n xnα, u3 (x) =

∞∑
n=0

u3,n xnα. (32)

Employing the initial conditions (30), then the kth approximate solution becomes as
follows:

u1k (x) = 1+

k∑
n=1

u1,n xnα, u2k (x) = 1+

k∑
n=1

u2,n xnα, u3k (x) = 1+

k∑
n=1

u3,n xnα. (33)



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 9 of 19

To get the kth approximate solution of System (29)-(30), we define the kth residual
functions of equations (29) as follows:

Rf (u1k (x)) = Dαu1k (x)− u1k (x) ,

Rf (u2k (x)) = Dαu2k (x)− 2u21k(x),

Rf (u3k (x)) = Dαu3k (x)− 3u1k (x)u2k (x) . (34)

After substituting the kth approximations (33) into Equation (34), then the kth residual
functions becomes as:

Rf (u1k (x)) =
k∑

n=1

u1,n
Γ (nα+ 1)

Γ ((n− 1)α+ 1)
x(n−1)α −

k∑
n=1

u1,n xnα − 1,

Rf (u2k (x)) =
k∑

n=1

u2,n
Γ (nα+ 1)

Γ ((n− 1)α+ 1)
x(n−1)α − 2

(
k∑

n=1

u1,n xnα + 1

)2

,

Rf (u3k (x)) =

k∑
n=1

u3,n
Γ (nα+ 1)

Γ ((n− 1)α+ 1)
x(n−1)α − 3

(
1 +

k∑
n=1

u1,n xnα

)

×

(
1 +

k∑
n=1

u2,n xnα

)
. (35)

According to the formulation used in the previous section, the kth approximations can
be achieved by obtaining the coefficients u1,j , u2,j , u3,j for j = 1, . . . , k via solving the
following algebraic equations iteratively:

lim
x→0

Rf
(
u1j (x)

)
x(j−1)α

= 0 , lim
x→0

Rf
(
u2j (x)

)
x(j−1)α

= 0, lim
x→0

Rf
(
u3j (x)

)
x(j−1)α

= 0 , j = 1, . . . , k.

(36)
These coefficients are summarized in Table 2.

So, the LRF solution of the System (29)-(30) has the following series expansions:

u1 (x) = 1 +
xα

Γ (α+ 1)
+

x2α

Γ (2α+ 1)
+

x3α

Γ (3α+ 1)
+

x4α

Γ (4α+ 1)
+ . . . ,

u2 (x) = 1 +
2xα

Γ (α+ 1)
+

4x2α

Γ (2α+ 1)
+

(
4

Γ (3α+ 1)
+

2 Γ (2α+ 1)

Γ 2 (α+ 1)Γ (3α+ 1)

)
x3α,

+

(
4

Γ (4α+ 1)
+

4 Γ (3α+ 1)

Γ (α+ 1)Γ (2α+ 1)Γ (4α+ 1)

)
x4α + . . . ,

u3 (x) = 1 +
3xα

Γ (α+ 1)
+

9x2α

Γ (2α+ 1)
+

(
15

Γ (3α+ 1)
+

6 Γ (2α+ 1)

Γ 2 (α+ 1)Γ (3α+ 1)

)
x3α



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 10 of 19

Table 2: The coefficients of the 4th approximations of u1 (x) , u2 (x) , u3 (x) for the System
(29)-(30).

k u1,k u2,k u3,k

0 1 1 1

1 1
Γ (α+1)

2
Γ (α+1)

3
Γ (α+1)

2 1
Γ (2α+1)

4
Γ (2α+1)

9
Γ (2α+1)

3 1
Γ (3α+1)

4
Γ (3α+1) +

2 Γ(2+1)
Γ2(α+1)Γ (3α+1)

15
Γ (3α+1) +

6 Γ(2+1)
Γ2(+1)Γ (3α+1)

4 1
Γ (4α+1)

4
Γ (4α+1) +

4 Γ(3+1)
Γ (α+1)Γ (2α+1)Γ (4α+1)

15
Γ (4α+1) +

6 Γ(2α+1)
Γ2(α+1)Γ (4α+1)

+ 18 Γ (3α+1)
Γ (α+1)Γ (2α+1)Γ (4α+1)

+

(
15

Γ (4α+ 1)
+

6 Γ (2α+1)

Γ 2 (α+ 1)Γ (4α+ 1)
+

18 Γ (3α+ 1)

Γ (α+ 1)Γ (2α+ 1)Γ (4α+ 1)

)
x4α + . . . .

(37)
For α = 1, the expansions in (37) become as follows:

u1 (x) = 1 + x+
x2

2
+

x3

3!
+

x4

4!
+ . . . ,

u2 (x) = 1 + 2x+ 2x2 +
4x3

3
+

2x4

3
+ . . . ,

u3 (x) = 1 + 3x+
9x2

2
+

9x3

2
+

27x4

8
+ . . . , (38)

which are the expansions of the exact solution indicated in (31).
Figure 2 shows the graphs of the 4th approximate solution of the system (29)-(30) at

different values of α, as well as the exact solution at α = 1. The graph shows agreement
between the approximate and actual solutions at α = 1. It also shows the effect of the
derivative’s order on the solution’s behavior. It is noticeable that the convergence intervals
of the solution decrease due to the action and effects of nonlinearity in the equation defining
the dependent variables.



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 11 of 19

(a) (b)

(c)

Figure 2: The curves of the 4th approximate and exact (solid curve) solutions of System
(29)-(30) at different values of α.

Example 4.3. Consider the following system of non-linear fractional differential equations

Dαu1 (x) = −1002 u1 (x) + 1000 u22 (x) ,

Dαu2 (x) = u1 (x)− u2 (x)− u22 (x) , (39)

where 0 < α ≤ 1, 0 ≤ x with the initial conditions

u1 (0) = 1, u2 (0) = 1. (40)

In the classical case (α = 1), the exact solution is given as follows:

u1 (x) = e−2x, u2 (x) = e−x. (41)

According to the methodology of the LRF method, we assume the solution of the
System (39)-(40) has fractional series expansions as:

u1 (x) =
∞∑
n=0

u1,n xnα, u2 (x) =
∞∑
n=0

u2,n xnα. (42)



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 12 of 19

Based on the initial conditions in Equation (40), the kth approximate solution can be
written as

u1k (x) = 1 +
k∑

n=1

u1,n xnα, u2k (x) = 1 +
k∑

n=1

u2,n xnα. (43)

The residual and kth residual functions of Equation (39) can be given respectively as
follows:

Rf (u1 (x)) = Dαu1 (x) + 1002 u1 (x)− 1000 u22 (x) ,

Rf (u2 (x)) = Dαu2 (x)− u1 (x) + u2 (x) + u22 (x) . (44)

Rf (u1k (x)) = Dαu1k (x) + 1002 u1k (x)− 1000 u22k(x),

Rf (u2k (x)) = Dαu2k (x)− u1k (x) + u2k (x) + u22k(x). (45)

Substituting the kth approximation (43) into the kth residual function (45) gives the series
form of the kth residual function as follows:

Rf (u1k (x)) =

k∑
n=1

u1,n
Γ (nα+ 1)

Γ ((n− 1)α+ 1)
x(n−1)α + 1002

(
1 +

k∑
n=1

u1,n xnα

)

−1000

(
1 +

k∑
n=1

u2,n xnα

)2

,

Rf (u2k (x)) =
k∑

n=1

u2,n
Γ (nα+ 1)

Γ ((n− 1)α+ 1)
x(n−1)α −

k∑
n=1

u1,n xnα +
k∑

n=1

u2,n xnα

+

(
1 +

k∑
n=1

u2,n xnα

)2

. (46)

According to the formulation used in the previous section, the kth approximations can be
achieved by obtaining the coefficients u1,j , u2,j for j = 1, . . . , k via solving the following
algebraic equations iteratively:

lim
x→0

Rf
(
u1j (x)

)
x(j−1)α

= 0 , lim
x→0

Rf
(
u2j (x)

)
x(j−1)α

= 0. (47)

These coefficients are summarized in Table 3.
So, the LRF solution of the System (39)-(40) has the following series expansions:

u1 (x) = 1− 2xα

Γ (α+ 1)
+

4x2α

Γ (2α+ 1)
+

(
− 2008

Γ (3α+ 1)
+

1000 Γ (2α+ 1)

Γ 2 (α+ 1)Γ (3α+ 1)

)
x3α

+

(
2014016

Γ (4α+ 1)
− 1004000 Γ (2α+ 1)

Γ 2 (α+ 1)Γ (4α+ 1)
− 2000 Γ (3α+ 1)

Γ (α+ 1)Γ (2α+ 1)Γ (4α+ 1)

)
x4α + . . . ,



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 13 of 19

Table 3: The coefficients of the 4th approximations of u1 (x) , u2 (x) for the System
(39)-(40).

k u1,k u2,k

0 1 1

1 − 2
Γ (α+1) − 1

Γ (α+1)

2 4
Γ (2α+1)

1
Γ (2α+1)

3 − 2008
Γ (3α+1) +

1000 Γ (2α+1)
Γ2(α+1)Γ (3α+1)

1
Γ (3α+1) −

Γ (2α+1)
Γ2(α+1)Γ (3α+1)

4 2014016
Γ (4α+1) −

1004000 Γ (2α+1)
Γ2(α+1)Γ (4α+1)

−2011
Γ (4α+1) +

1003 Γ (2α+1)
Γ2(α+1)Γ (4α+1)

− 2000 Γ (3α+1)
Γ (α+1)Γ (2α+1)Γ (4α+1) + 2 Γ (3α+1)

Γ (α+1)Γ (2α+1)Γ (4α+1)

u2 (x) = 1− xα

Γ (α+ 1)
+

x2α

Γ (2α+ 1)
+

(
1

Γ (3α+ 1)
− Γ (2α+ 1)

Γ 2 (α+ 1)Γ (3α+ 1)

)
x3α

+

(
− 2011

Γ (4α+ 1)
+

1003 Γ (2α+ 1)

Γ 2 (α+ 1)Γ (4α+ 1)
+

2 Γ (3α+ 1)

Γ (α+ 1)Γ (2α+ 1)Γ (4α+ 1)

)
x4α + . . . .

(48)

(a) (b)

Figure 3: The curves of the 4th approximate and exact (solid curve) solutions of System
(39)-(40) at different values of α.

For α = 1, the expansions in (48) become as follows:

u1 (x) = 1− 2x+ 2x2 − 4x3

3
+

2x4

3
+ . . . ,

u2 (x) = 1− x+
x2

2
− x3

6
+

x4

24
+ . . . . (49)

which are the expansions of the exact solution indicated in (41).
Figure 3 shows graphs of the 4th approximate solution of the system (39)-(40) for

different values of α. The graphs show the intervals of convergence and the effect of the



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 14 of 19

order of the derivative on the behavior of the solution.

Example 4.4. Consider the following system of non-linear fractional differential equations

Dαu1 (x) = u21 (x) + u2 (x) ,

Dαu2 (x) = u2 (x) cos (u1 (x)) , (50)

where 0 < α ≤ 1, x ≥ 0, with the initial conditions

u1 (0) = 0, u2 (0) = 1. (51)

According to the methodology of the LRF method, we assume the solution of the
System (50)-(51) has fractional PS expansions as:

u1 (x) =
∞∑
n=0

u1,n xnα, u2 (x) =
∞∑
n=0

u2,n xnα. (52)

Based on the initial conditions in Equation (51), the kth approximate solution can be
written as

u1k (x) =

k∑
n=1

u1,n xnα, u2k (x) = 1 +

k∑
n=1

u2,n xnα. (53)

The residual and kth residual functions of Equation (50) can be given respectively as
follows:

Rf (u1 (x)) = Dαu1 (x)− u21 (x)− u2 (x) ,

Rf (u2 (x)) = Dαu2 (x)− u2 (x) cos (u1 (x)) . (54)

Rf (u1k (x)) = Dαu1k (x)− u21k (x)− u2k (x) ,

Rf (u2k (x)) = Dαu2k (x)− u2k (x) cos (u1k (x)) . (55)

By substituting the first approximations u11 (x) = u1,1 xα, u21 (x) = 1 + u2,1x
α into the

first residual functions, Rf (u11 (x)) , Rf (u21 (x)) , we obtain

Rf (u11 (x)) = u1,1Γ (α+ 1)− (u1,1x
α)2 − 1− u2,1x

α,

Rf (u21 (x)) = u2,1Γ (α+ 1)− (1 + u2,1x
α) cos (u1,1 xα) . (56)

Solving the equations limx→0Rf (u11 (x)) = 0, limx→0Rf (u21 (x)) = 0 for u1,1, u2,1,
respectively, we have

u1,1 =
1

Γ (α+ 1)
, u2,1 =

1

Γ (α+ 1)
. (57)

In order to determine the second coefficients u1,2 and u2,2, substitute the second ap-

proximations u12 (x) =
xα

Γ (α+1) +u1,2x
2α, u22 (x) = 1+ xα

Γ (α+1) +u2,2x
2α into Rf (u12 (x)) ,

Rf (u22 (x)) to have

Rf (u12 (x)) = 1+ u1,2
Γ (2α+ 1) xα

Γ (α+ 1)
−
(

xα

Γ (α+ 1)
+ u1,2 x2α

)2

−
(
1 +

xα

Γ (α+ 1)
+ u2,2 x2α

)
,



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 15 of 19

Rf (u22 (x)) = 1 + u2,2
Γ (2α+ 1) xα

Γ (α+ 1)
−
(
1 +

xα

Γ (α+ 1)
+ u2,2x

2α

)
cos

(
xα

Γ (α+ 1)
+ u1,2x

2α

)
.

(58)

Making simple calculations, the limits limx→0
Rf(u12(x))

xα = 0, limx→0
Rf(u22(x))

xα = 0, yield

u1,2 =
1

Γ (2α+ 1)
, u2,2 =

1

Γ (2α+ 1)
. (59)

The values of the coefficients u1,3 and u2,3 are also obtained by substituting u13 (x) =
xα

Γ(α+1) +
2

Γ(2α+1)x
2α + u1,3 x3α and u23 (x) = 1 + xα

Γ(α+1) + u2,3 x3α into Rf (u13 (x)) , Rf (u23 (x)), and

then we solve the following equations:

lim
x→0

Rf (u13 (x))

x2α
= 0, lim

x→0

Rf (u23 (x))

x2α
= 0. (60)

Then we have

u1,3 =
Γ 2 (α+ 1) + Γ (2α+ 1)

Γ 2 (α+ 1)Γ (3α+ 1)
, u2,3 =

2Γ 2 (α+ 1)− Γ (2α+ 1)

2Γ 2 (α+ 1)Γ (3α+ 1)
. (61)

Similarly, we can find the following coefficients, which are as follows:

u1,4 =
2Γ 2 (α+ 1)Γ (2α+ 1)− Γ 2 (2α+ 1) + 4Γ (α+ 1)Γ (3α+ 1)

2Γ 2 (α+ 1)Γ (2α+ 1)Γ (4α+ 1)
,

u2,4 =
Γ (α+ 1) (2Γ 2 (α+ 1)− Γ (2α+ 1))Γ (2α+ 1)− (2Γ 2 (α+ 1) + Γ (2α+ 1))Γ (3α+ 1)

2Γ 3 (α+ 1)Γ (2α+ 1)Γ (4α+ 1)
.

(62)
So, the LRF solution of the System (50)-(51) has the following series expansions:

u1 (x) =
xα

Γ (α+ 1)
+

x2α

Γ (2α+ 1)
+
Γ 2 (α+ 1) + Γ (2α+ 1)

Γ 2 (α+ 1)Γ (3α+ 1)
x3α

+
2Γ 2 (α+ 1)Γ (2α+ 1)− Γ 2 (2α+ 1) + 4Γ (α+ 1)Γ (3α+ 1)

2Γ 2 (α+ 1)Γ (2α+ 1)Γ (4α+ 1)
x4α + . . . ,

u2 (x) = 1+
xα

Γ (α+ 1)
+

x2α

Γ (2α+ 1)
+
2Γ 2 (α+ 1)− Γ (2α+ 1)

2Γ 2 (α+ 1)Γ (3α+ 1)
x3α

+
Γ (α+ 1) (2Γ 2 (α+ 1)− Γ (2α+ 1))Γ (2α+ 1)− (2Γ 2 (α+ 1) + Γ (2α+ 1))Γ (3α+ 1)

2Γ 3 (α+ 1)Γ (2α+ 1)Γ (4α+ 1)
x4α+. . . .

(63)
For α = 1, the solution (63) becomes as follows:

u1 (x) = x+
x2

2
+

x3

2
+

x4

4
+ . . . ,

u2 (x) = 1 + x+
x2

2
− x4

4
+ . . . . (64)

This solution coincides with the solution attained by using the Adomian decomposition method
[38].



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 16 of 19

Tables 4 and 5 present numerical values for the 4th approximation of the solution to the IVP
(50)-(51) and the residual error of the approximate solution at various α values. Since there is no
exact solution available for this problem, we focus solely on the residual error, which is defined for
the problem as follows:

Rf. Err. (x) =
∣∣Dαu1k (x)− u2

1k
(x)− u2k (x)

∣∣ ,
Rf. Err. (x) = |Dαu2k (x)− u2k (x) cos (u1k (x)) | . (65)

Table 4: The 4th approximation of u1 (x), the IVP (33)-(34) solution, and the residual
error for α = 1 and α = 0.8.

α = 1 α = 0.8

t u14 (x) Rf. Err. u14 (x) Rf. Err.

0.0 0 0 0 0

0.1 0.105525 1.10526×10−4 0.191650 1.47602×10−3

0.2 0.224400 1.95536×10−3 0.371681 1.58855×10−2

0.3 0.360525 1.09533×10−2 0.573817 6.71247×10−2

0.4 0.518400 3.83386×10−2 0.807786 1.92998×10−1

0.5 0.703125 1.0376×10−1 1.080985 4.48851×10−1

The residual error serves as an indicator of the solution’s accuracy, although it does not precisely
measure the error in the same way that absolute error does. Even though the residual error remains
within a range of 1 to 3 and is concentrated in a small interval, the solution period can be extended
by employing alternative techniques, such as multi-step methods.

Table 5: The 4th approximation of u2 (x), the IVP (33)-(34) solution, and the residual
error for α = 1 and α = 0.8.

α = 1 α = 0.8

t u24 (x) Rf. Err. u24 (x) Rf. Err.

0.0 1 0 1 0
0.1 1.104975 1.71532×10−4 1.187653 2.13071×10−3

0.2 1.219600 2.97806×10−3 1.347857 2.22820×10−2

0.3 1.342975 1.63625×10−2 1.504286 9.14918×10−2

0.4 1.473600 5.60118×10−2 1.657010 2.53724×10−1

0.5 1.609375 1.47328×10−1 1.803767 5.60315×10−1

5. Conclusion

In this article, we introduced a new analytical iterative technique called the LRF method,
which we used to solve linear and nonlinear systems of Caputo-fractional differential equations.
This approach simplifies discovering exact solution patterns while reducing the need for complex
computational calculations. The main advantage of this method is its straightforwardness in
computing the coefficients of the series solution by evaluating the limit of a form that includes



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 17 of 19

the residual function for the given equations. This is a departure from other well-known analytic
techniques that rely on differential and integral operators, which can be challenging in the fractional
case. We plan to modify our approach to handle more complex real-world applications in future
work.

Acknowledgements

The authors express their gratitude to the dear referees, who wish to remain anonymous, and
the editor for their helpful suggestions that improved the final version of this paper. This research
is funded by Zarqa University-Jordan.

References

[1] A. Afreen and A. Raheem. Study of a nonlinear system of fractional differential equations with
deviated arguments via adomian decomposition method. International Journal of Applied and
Computational Mathematics, 8(5):269, 2022.

[2] K. I. Ahmed, H. D. Adam, N. Almutairi, and S. Saber. Analytical solutions for a class of
variable-order fractional liu system under time-dependent variable coefficients. Results in
Physics, 56:107311, 2024.

[3] M. Alaroud. Application of laplace residual power series method for approximate solutions of
fractional ivp’s. Alexandria Engineering Journal, 61(2):1585–1595, 2022.

[4] N. Almutairi and S. Saber. On chaos control of nonlinear fractional newton-leipnik system
via fractional caputo-fabrizio derivatives. Scientific Reports, 13(1):22726, 2023.

[5] N. Almutairi and S. Saber. Application of a time-fractal fractional derivative with a power-
law kernel to the burke-shaw system based on newton’s interpolation polynomials. MethodsX,
12:102510, 2024.

[6] M. Altalla, B. Shanmukha, A. El-Ajou, and M. Alkord. Taylor’s series in terms of the mod-
ified conformable fractional derivative with applications. Nonlinear Functional Analysis and
Applications, 2024:435–450, 2024.

[7] N. Anjum, C. H. He, and J. H. He. Two-scale fractal theory for the population dynamics.
Fractals, 29(07):2150182, 2021.

[8] A. Atangana and D. Baleanu. New fractional derivatives with nonlocal and non-singular
kernel: theory and application to heat transfer model. Thermal Science, 20:763–769, 2016.

[9] A. Burqan. A novel scheme of the ara transform for solving systems of partial fractional
differential equations. Fractal and Fractional, 7(4):306, 2023.

[10] A. Burqan, M. Khandaqji, Z. Al-Zhour, A. El-Ajou, and T. Alrahamneh. Analytical approx-
imate solutions of caputo fractional kdv-burgers equations using laplace residual power series
technique. Journal of Applied Mathematics, 1:7835548, 2024.

[11] A. Burqan, M. Shqair, A. El-Ajou, S. Ismaeel, and Z. Al-Zhour. Analytical solutions to the
coupled fractional neutron diffusion equations with delayed neutrons system using laplace
transform method. AIMS Mathematics, 8(8):19297–19312, 2023.

[12] M. Caputo. Linear models of dissipation whose q is almost frequency independent–ii. Geo-
physical Journal International, 13(5):529–539, 1967.

[13] S. Cifani and E. R. Jakobsen. Entropy solution theory for fractional degenerate convection–
diffusion equations. Annales de l’IHP Analyse non linéaire, 28(3):413–441, 2011.

[14] S. Das. Analytical solution of a fractional diffusion equation by variational iteration method.
Computers Math. Appl., 57(3):483–487, 2009.

[15] A. Dawar, H. Khan, S. Islam, and W. Khan. The improved residual power series method for



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 18 of 19

a system of differential equations: a new semi-numerical method. International Journal of
Modelling and Simulation, pages 1–14, 2023.

[16] K. Diethelm, N. J. Ford, and A. D. Freed. A predictor-corrector approach for the numerical
solution of fractional differential equations. Nonlinear Dynamics, 29:3–22, 2002.

[17] A. El-Ajou. Taylor’s expansion for fractional matrix functions: theory and applications.
Mathematics, 21(7):1–17, 2020.

[18] A. El-Ajou and A. Burqan. Limit residual function method and applications to pde models.
The European Physical Journal Plus, 139(11):973, 2024.

[19] A. El-Ajou and A. Burqan. A new algorithm for generating power series solutions for a broad
class of fractional pdes: Applications to interesting problems. Fractals, 2450120, 2024.

[20] A. El-Ajou, M. Shqair, I. Ghabar, A. Burqan, and R. Saadeh. A solution for the neutron
diffusion equation in the spherical and hemispherical reactors using the residual power series.
Frontiers in Physics, 11:1229142, 2023.

[21] J. He. Homotopy perturbation method: a new nonlinear analytical technique. Applied Math-
ematics and Computation, 135(1):73–79, 2003.

[22] J. He, M. Jiao, K. A. Gepreel, and Y. Khan. Homotopy analysis method: A new analytical
technique for nonlinear problems. Mathematics and Computers in Simulation, 204:243–258,
2023.

[23] J. He and X. Wu. Variational iteration method: New development and applications. Com-
puters & Mathematics with Applications, 54(7-8):881–894, 2007.

[24] J. H. He, Z. B. Li, and Q. L. Wang. A new fractional derivative and its application to
explanation of polar bear hairs. Journal of King Saud University-Science, 28(2):190–192,
2016.

[25] J. H. He, L. Verma, B. Pandit, A. K. Verma, and R. P. Agarwal. A new taylor series based
numerical method: Simple, reliable, and promising. Journal of Applied and Computational
Mechanics, 9(4):1122–1134, 2023.

[26] H. Jafari and V. Daftardar-Gejji. Solving a system of nonlinear fractional differential equa-
tions using adomian decomposition. Journal of Computational and Applied Mathematics,
196(2):644–651, 2006.

[27] H. Khresat, A. El-Ajou, S. Al-Omari, S. E. Alhazmi, and M. N. Oqielat. Exact and approxi-
mate solutions for linear and nonlinear partial differential equations via laplace residual power
series method. Axioms, 12(7):694, 2023.

[28] W. Li and Y. Pang. Application of adomian decomposition method to nonlinear systems. Adv
Differ Equ, 67:2020, 2020.

[29] R. L. Magin, C. Ingo, L. Colon-Perez, W. Triplett, and T. H. Mareci. Characterization of
anomalous diffusion in porous biological tissues using fractional order derivatives and entropy.
Microporous and Mesoporous Materials, 178:39–43, 2013.

[30] B. A. Mahmood, M. A. Yousif, and L. Liu. A residual power series technique for solving
boussinesq–burgers equations. Cogent Mathematics, 4(1), 2017.

[31] F. Mainardi. Fractional calculus: Some basic problems in continuum and statistical mechanics.
In A. Carpinteri and F. Mainardi, editors, Fractals and Fractional Calculus in Continuum
Mechanics, pages 291–348. Springer-Verlag, 1997.

[32] F. Mainardi, M. Raberto, R. Gorenflo, and E. Scalas. Fractional calculus and continuous-time
finance ii: the waiting-time distribution. Physica A: Statistical Mechanics and its Applications,
287(3-4):468–481, 2000.

[33] R. K. Nagle, E. B. Saff, and A. D. Snider. Fundamentals of differential equations and boundary
value problems. Pearson Custom Publishing, Boston, Mass., 2012.

[34] A. Saadatmandi and M. Dehghan. A new operational matrix for solving fractional-order



A. El-Ajou, A. Burqan / Eur. J. Pure Appl. Math, 18 (1) (2025), 5617 19 of 19

differential equations. Computers Math. Appl., 59(3):1326–1336, 2010.
[35] S. Saber. Control of chaos in the burke-shaw system of fractal-fractional order in the sense of

caputo-fabrizio. Journal of Applied Mathematics and Computational Mechanics, 23(1):83–96,
2024.

[36] A. Sarhan, A. Burqan, R. Saadeh, and Z. Al-Zhour. Analytical solutions of the nonlinear time-
fractional coupled boussinesq-burger equations using laplace residual power series technique.
Fractal and Fractional, 6(11):631, 2022.

[37] S. Zhang and H. Q. Zhang. Fractional sub-equation method and its applications to nonlinear
fractional pdes. Physics Letters A, 375(7):1069–1073, 2011.

[38] E. A. A. Ziada. Analytical solution of nonlinear system of fractional differential equations.
Journal of Applied Mathematics and Physics, 9(10):2544–2557, 2021.