EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 17, No. 4, 2024, 2939-2961
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Multi-step Residual Power Series Method for Solving
Stiff Systems

Mohammad Al Zurayqat1, Shatha Hasan1,2,∗

1 Department of Applied Science, Ajloun College, Al-Balqa Applied University,
Ajloun 26816, Jordan
2 Jadara Research Center, Jadara University, Irbid 21110, Jordan

Abstract. In this paper, an efficient algorithm based on the residual power series method (RPSM)
is presented to solve stiff systems of Caputo fractional order. We apply the RPSM on subintervals
to get approximate solutions of these types of systems. The RPSM has advantages that it is suitable
to solve linear and nonlinear systems and it gives high accurate results. Modifying this technique
to multi-step RPSM considerably reduces the number of arithmetic operations and so reduces the
time, especially when dealing with Stiff systems. Several numerical examples are given to show the
efficiency, simplicity and the accuracy of the proposed method. Comparing classical RPSM with
the new multi-step scheme shows that multi-step RPSM controls the convergence behaviour of the
stiff systems. That is, the comparison reveals that MS-FRPSM reduces both absolute and residual
errors. More iterations and a smaller step size lead to higher accuracy. Moreover, in MS-FRPSM,
the intervals of convergence for the series solution will increase.

2020 Mathematics Subject Classifications: 26A33, 34A08, 65L04

Key Words and Phrases: Multi-step method, Residual power series method, Stiff system,
Numerical solution, Caputo fractional derivative

1. Introduction

In various fields such as chemical kinetics, aerodynamics, ballistics, electrical circuit
theory, and missile guidance, there are specific types of differential equations that pose
significant challenges for solution using classical numerical methods [6]. As a result, many
analytic and numerical methods have been developed to deal with such equations. Exam-
ples for these methods and the real-life applications of differential equations can be found
in [14], [15],[23],[1], and [2].
Curtiss and Hirschfelder were the first to bring attention to these challenging differential
equations, commonly referred to as stiff equations [8]. The mathematical stiffness of a
problem arises from the disparate rates of various processes within the considered physical

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v17i4.5437

Email addresses: Mohammadzraqat22@gmail.com (M. Al Zurayqat), shatha@bau.edu.jo (S. Hasan)

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



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2940

models. This phenomenon occurs when certain solution components decay significantly
faster than others, characterized by terms such as e−λt where λ > 0.
In the last few years, several analytical methods have been exploited to give approximate
solutions for stiff systems of ordinary differential equations such as homotopy perturba-
tion method [4], block method [3], second derivative multistep method [22], variational
iteration method [5], and multi-step reproducing kernel Hilbert space method [12].
The residual power series method is one of the efficient numerical techniques employed
for solving various types of differential equations. For a comprehensive understanding
of this method and its applications in solving differential equations as well as integro-
differential equations of distinct nature, readers are encouraged to refer to the relevant
literature[16],[11], [10],[13], [21],[17], [20].

Unfortunately, when attempting to employ this method to tackle a stiff system, it
encounters limitations. The resulting approximate solutions remain valid only within a
narrow interval, displaying slow convergence or complete divergence when applied to a
broader range. It gives numerical results with rapidly increasing error of approximation.
Moreover, increasing the number of iterations results in a necessity of a large computer
memory and large time to run. To overcome these drawbacks, we formulate an algorithm
that depends on dividing any time interval into small subintervals and applying the RPSM
on each subinterval. We call this technique a multi-step residual power series method
(MS- RPSM). However, fractional differential equations (FDEs) become widely referred
because of their capabilities to model several problems in engineering and science such
as mechanical systems, dynamical systems, heat transfer, control theory, mixed convec-
tion flows, unification of diffusion, image processing, and wave propagation phenomenon
[18],[7],[19],[25], and [24]. This is because the fascinating properties of distinct fractional
operators, such as the ability to capture long-range dependencies and memory effects in
various systems. As a result, we apply in this work the MS-RPSM as a treatment to get
accurate solutions for stiff systems with the well-known Caputo fractional derivative.

2. Basics of Caputo operators and RPSM

Fractional residual power series method (FRPSM) is an analytical as well as numerical
method based on the concept on the error function and the generalized Taylor series. In
fact, it produces solutions in the form of rapidly convergent fractional power series (FPS)
with computationally straightforward components. Basic definitions that are related to
FPS and Caputo derivative are given below.

Definition 1. [9] A power series representation of the form

∞∑
k=0

bk (z − z0)
kβ = b0 + b1 (z − z0)

β + b2 (z − z0)
2β + · · · ,

where 0 ≤ m− 1 < β ≤ m ∈ N, is called FPS about z = z0 and bk, k = 0, 1, 2, . . . are the
coefficients for the series.
In order to present the definition of Caputo derivative, we need to define the space Cn

γ [t0, b]



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2941

as follows: Let q(t) be a function defined on [t0, b] and γ be a real number. We say that
q(t) is in the space Cγ [t0, b] if there is a continuous function q1(t) on [t0, b] such that
q(t) = tρq1(t) for a real number ρ > Cγ . Moreover, it is said to be in the space Cn

γ [t0, b] if
its classical nth derivative is in the space Cγ [t0, b].

Definition 2. [19] Assume that n − 1 < β ≤ n, n ∈ N, and q(t) ∈ Cn
−1[t0, b]. Then the

Caputo derivative of q(t) of order β is given by

C
t0D

β
t q(t) =

{
1

Γ(n−β)

∫ t
t0

dnq(τ)
dτn (t− τ)n−β−1 dτ, t > t0,

dnq(t)
dtn , β = n.

(1)

Now, we summarize the FRPSM as described in [10] through the following algorithm.
Algorithm 2.1 To get the FPS solution of the fractional initial value problem (FIVP)

C
t0D

β
t q(t) = B(t, q(t)), t ≥ t0, 0 < β ≤ 1, (2)

q(t0) = q0,

where B is a linear or nonlinear function of t and the unknown q(t), and q0 ∈ R is the
initial condition, we follow the steps:
Step 1: Assume the solution of (2) is of the FPS form

q(t) =

∞∑
v=0

bv(t− t0)
vβ. (3)

Step 2: Use the initial condition q(t0) = q0 as the zeroth coefficient of the FPS. That is,
b0 = q0.
Step 3: Define the R-th truncated FPS as

qR(t) = b0 +
R∑

v=1

bv(t− t0)
vβ. (4)

Step 4: Define the residual function and the N -th residual function, respectively, as

Residq(t) =
C
0 D

β
t q(t)−B(t, q(t)), (5)

Residq,N (t) = C
0 D

β
t qN (t)−B(t, qN (t)). (6)

Obviously,
Residq(t) = lim

N→∞
Residq,N (t) = 0,

and
C
0 D

β
t Residq(t) = 0, for t ≥ 0.

Step 5: For k = 0, 1, 2, . . . , R− 1, compute

(C0 D
(k−1)β
t Residq,N )(t0) = 0.

Step 6: Solve the resulting equation in Step 5 to get the coefficient bR.
Step 7: Repeat Steps 3-6 N -times to get the N -th truncated FPS where N can be chosen
to achieve the required accuracy.



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2942

3. Description of the MS-RPSM

In this section, we provide a modification to Algorithm 2.1 to get the MS-FRPSM in
order to solve stiff systems of Caputo fractional order, which may be of the specific form

C
t D

β
aGi(t) = Fi(t, G1, G2, . . . , Gm), Gi(a) = ci, (7)

a ≥ 0, t ∈ [a, b], i = 1, 2, . . . ,m, m ∈ N, 0 < β ≤ 1.
To clarify the MS-RPSM in which we will use to obtain the approximate solution of
system (7), we split the interval [a,b] into subintervals [tn−1, tn],n = 1, 2, . . . ,M , of step
size h = b−a

M and nodes t0 = a, tn = a+nh. Then we solve each FIVP on its corresponding
subinterval:

C
t D

β
tn−1

nGi(t) = Fi (t,
nG1,

nG2, . . . ,
nGm) , 0 < β ≤ 1, (8)

1Gi(t0) = ci,
nGi(tn−1) =

n−1Gi(tn−1),

t ∈ [tn−1, tn], i = 1, 2, . . . ,m, m ∈ N.

Now, we assume the solution of system (7) has the form:

Gi(t) =


1Gi(t), if t ∈ [a, t1],
2Gi(t), if t ∈ [t1, t2],
...

MGi(t), if t ∈ [tM−1, b],

=

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

∑∞
v=0

av1(t−t0)vβ

Γ(vβ+1) if t ∈ [a, t1],∑∞
v=0

av2(t−t1)vβ

Γ(vβ+1) if t ∈ [t1, t2],
...∑∞
v=0

avM (t−tM−1)
vβ

Γ(vβ+1) if t ∈ [tM−1, b],

and the N-th truncated solution has the form

Gi,N (t) =


1Gi,N (t), if t ∈ [a, t1],
2Gi,N (t), if t ∈ [t1, t2],
...

MGi,N (t), if t ∈ [tM−1, b],

=

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

∑N
v=0

av1(t−t0)vβ

Γ(vβ+1) if t ∈ [a, t1],∑N
v=0

av2(t−t1)vβ

Γ(vβ+1) if t ∈ [t1, t2],
...∑N
v=0

avM (t−tM−1)
vβ

Γ(vβ+1) if t ∈ [tM−1, b],

for each i = 1, 2, . . . ,m. Hence the residual function can be written as follows:

ResidGi(t) =


C
t D

β
a
1Gi(t)− Fi(t,

1G1,
1G2, . . . ,

1Gm), if t ∈ [a, t1],
C
t D

β
t1

2Gi(t)− Fi(t,
2G1,

2G2, . . . ,
2Gm), if t ∈ [t1, t2],

...
C
t Dt

β
M−1

MGi(t)− Fi(t,
MG1,N ,MG2,N , . . . ,MGm,N ), if t ∈ [tM−1, b],



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2943

and the Nth residual function has the form:

ResidGi,N
(t) =


C
t D

β
a
1Gi,N (t)− Fi(t,

1G1,N ,1G2,N , . . . ,1Gm,N ), t ∈ [a, t1]
C
t D

β
t1

2Gi,N (t)− Fi(t,
2G1,N ,2G2,N , . . . ,2Gm,N ), t ∈ [t1, t2]

...
...

C
t D

β
tM−1

MGi,N (t)− Fi(t,
M G1,N ,M G2,N , . . . ,M Gm,N ), t ∈ [tM−1, b]

for each i = 1, 2, . . . ,m.
For more illustrations, we apply the FRPSM to the system with n = 1. That is, to

the first subinterval on which the system has the form:

C
t D

β

a
1Gi(t) = Fi(t,

1G1,
1G2, . . . ,

1Gm), 0 < β ≤ 1, t ∈ [a, t1],
1Gi(a) = ci, (9)

to get the Nth truncated FPS of the form:

1Gi,N (t) =

N∑
v=0

av1(t− a)vβ

Γ(vβ + 1)
, t ∈ [a, t1], i = 1, 2, . . . ,m.

After that, we take the second system, i.e., when n = 1, which has the form

C
t D

β

t1
2Gi(t) = Fi(t,

2G1,
2G2, . . . ,

2Gm), 0 < β ≤ 1, t ∈ [t1, t2],
2Gi(t1) =

1 Gi(t1).
(10)

Applying the FRPSM on system (10) gives us the Nth truncated FPS:

2Gi,N (t) =

N∑
v=0

av2(t− t1)
vβ

Γ(vβ + 1)
, t ∈ [t1, t2], i = 1, 2, . . . ,m.

Repeating this process for the M systems, we get the approximate solution for the stiff
system in (7) in the form of a piecewise-defined function in which its components are the
Nth truncated FPS:

Gi,N (t) =


1Gi,N (t), t ∈ [a, t1],
2Gi,N (t), t ∈ [t1, t2],
...

...
MGi,N (t), t ∈ [tM−1, b],

and the exact solution can be obtained from limN→∞Gi,N (t).

4. Numerical applications

In this section, we test the efficiency of the MS-FRPSM in solving stiff system through
solving three numerical examples. We make a comparison between classical FRPS and
MS-FRPSM in some examples to see the importance of minimizing the interval length



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2944

to obtain better approximations. Since the exact solutions of fractional stiff systems are
not available, we compute the residual errors to check our computations. In all of our
examples, we used Mathematica 10 software.
Example 4.1 Consider the following linear fractional system:

C
t D

β
0G(t) = −G(t) + 95W (t), t ∈ [0, 1], 0 < β ≤ 1, (11)

C
t D

β
0W (t) = −G(t) + 97W (t),

G(0) = 1, W (0) = 1.

The exact solution for the system (11) when β = 1 is

G(t) =
1

47
e−96t(−48 + 95e94t), W (t) =

1

47
e−96t(−48 + e94t).

If we try to solve (11) using the classical FRPSM, we assume the solution has the FPS:

G(t) = 1 +
∞∑
v=1

av
Γ(1 + vβ)

tvβ, and W (t) = 1 +
∞∑
v=1

bv
Γ(1 + vβ)

tvβ.

And now we define the N -th truncated series for the forms

GN (t) = 1 +
N∑
v=1

av
Γ(1 + vβ)

tvβ, and WN (t) = 1 +
N∑
v=1

bv
Γ(1 + vβ)

tvβ.

The N -th truncated residual function are:

ResidG,N =C
t Dβ

0GN (t)+GN (t)−95WN (t), t ∈ [0, 1]ResidW,N =C
t Dβ

0WN (t)+GN (t)−97WN (t).

When N = 1, we solve ResidW,1(0) = 0 and ResidG,1(0) = 0 to get the values of a1 and b1
as:

a1 = 94, b1 = −98.

When N = 2, we solve C
0 D

β
t ResidW,2(0) = 0 and C

t D
β
0 ResidG,2(0) = 0 to get the values of

a2 and b2 as
a2 = −9404, b2 = 9412.

When N = 3, we solve C
t D

2β
0 ResidW,3(0) = 0 and C

t D
2β
0 ResidG,3(0) = 0 to get the values

of a3 and b3 as
a3 = 903544, b3 = −903560,

and so on up to theNth coefficients which can be obtained by solving C
t D

(N−1)β
0 ResidW,N (0) =

0 and C
t D

(N−1)β
0 ResidG,N (0) = 0. Some of the first coefficients are obtained as follows:

a4 = −86741744, b4 = 86741776,



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2945

a5 = 832710464, b5 = −832710528,

a6 = −799412210624, b6 = 79912210752,

a7 = 76743572232064, b7 = −76743572232320,

a8 = −7367382934302464, b8 = 7367382934302976,

a9 = 707268761693085184, b9 = −707268761693086208,

a10 = −67897801122536274944, b10 = 67897801122536276992.

In fact, we use Mathematica software to get the 30th truncated FPS of G(t) and W (t),
and compute some numerical and display some graphical results. Table 1 shows the val-
ues of the exact and approximate solutions when β = 1 in addition to the absolute error.
Also, it presents numerical solutions for different values of β. The residual error function
is computed in Table 2 for β = 1, 0.99, 0.98, 0.97 and the approximate solutions for the
same values of β are shown in Figure 1.

Table 1 Numerical results for classical FRPSM for example 4.1

ti G(t) G30(t), β = 1 |G(t)−G30(t)|, β = 1 G30(t), β = 0.99 G30(t), β = 0.98 G30(t), β = 0.97

0.0 1.00000 1.00000 0.00000 1.00000 1.00000 1.00000
0.1 1.65481 1.65454 0.000269072 1.64314 1.62541 1.57125
0.2 1.35490 −468189 468191 −2.19151× 106 −1.02179× 1024 −4.74525× 107

0.3 1.10930 −1.13099× 1011 1.13099× 1011 −4.64964× 1011 −1.90422× 1024 −7.76854× 1012

0.4 0.908218 −7.27568× 1014 7.27568× 1014 −2.72922× 1015 −1.01996× 1024 −3.79751× 1016

0.5 0.74359 −6.45204× 1017 6.45204× 1017 −2.25504× 1018 −7.85291× 1024 −2.72471× 1019

0.6 0.608797 −1.63816× 1020 1.63816× 1020 −5.40562× 1020 −1.77741× 1024 −5.82343× 1021

0.7 0.498441 −1.75739× 1022 1.75739× 1022 −5.52509× 1022 −1.73098× 1024 −5.40409× 1023

0.8 0.408089 −1.00451× 1024 1.004514× 1024 −3.02892× 1024 −9.10183× 1024 −2.72564× 1025

0.9 0.334115 −3.55209× 1025 3.55209× 1025 −1.03246× 1026 −2.99082× 1024 −8.63426× 1026

1.0 0.27355 −8.60381× 1026 8.60381× 1026 −2.42025× 1028 −6.78535× 1024 −1.89592× 1028

ti W (t) W30(t), β = 1 |W (t)−W30(t)|, β = 1 W30(t), β = 0.99 W30(t), β = 0.98 W30(t), β = 0.97

0.0 1.00000 1.00000 0.00000 1.00000 1.00000 1.00000
0.1 0.01735 - 0.0170816 0.000269072 −0.0143377 −0.00539508 0.0398389
0.2 0.0142621 468191 468191 2.19151× 106 1.02179× 107 4.74525× 107

0.3 0.0116768 1.13099× 1011 1.13099× 1011 4.64964× 1011 1.90422× 1012 7.76854× 1012

0.4 0.00956019 7.27568× 1014 7.27568× 1014 2.72922× 1015 1.01996× 1016 3.79751× 1016

0.5 0.00782722 6.45204× 1017 6.45204× 1017 2.25504× 1018 7.85291× 1018 2.72471× 1019

0.6 0.00640839 1.63816× 1020 1.63816× 1020 5.40562× 1020 1.77741× 1021 5.82343× 1021

0.7 0.00524674 1.75739× 1022 1.75739× 1022 5.52509× 1022 1.73098× 1023 5.40409× 1023

0.8 0.003517 1.00451× 1024 1.00451× 1024 3.02892× 1024 9.10183× 1024 −2.72564× 1025

0.9 0.003517 3.55209× 1025 3.55209× 1025 1.03246× 1026 2.99082× 1026 8.63426× 1026

1.0 0.00287947 8.60381× 1026 8.60381× 1026 2.42025× 1028 6.78535× 1024 1.89592× 1028

We note from Tables 1-2 and Figure 1 that the absolute and the residual errors are
very large and still rapidly increase. Hence, we have to apply the MS-FRPSM so we may
reduce the errors. We solve system (11) again with a step size of h = 0.1, so we divide the
interval into 10 subintervals. Assuming the solution has the form:



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2946

Table 2 Some numerical values of the residual function for classical FRPSM, example 4.1

ti |Resid(G,30)(t)| |Resid(W,30)(t)|
β = 1 β = 0.99 β = 0.98 β = 0.97 β = 1 β = 0.99 β = 0.98 β = 0.97

0.1 3.27× 10−3 3.00522 2.13794 2.21345 4.49× 10−3 2.10× 108 9.80× 108 4.55× 109

0.2 4.49× 107 2.10× 108 9.80× 108 4.55× 109 1.08× 1013 4.46× 1013 1.82× 1014 7.45× 1014

0.3 1.08× 1013 4.46× 1013 1.82× 1014 7.45× 1014 6.98× 1016 2.62× 1017 9.79× 1017 3.64× 1018

0.4 6.98× 1016 2.62× 1017 9.79× 1017 3.64× 1018 6.19× 1019 2.16× 1020 7.53× 1020 2.61× 1021

0.5 6.19× 1019 2.16× 1020 7.53× 1020 2.61× 1021 1.57× 1022 5.18× 1022 1.70× 1023 5.59× 1023

0.6 1.57× 1022 5.18× 1022 1.70× 1023 5.59× 1023 1.68× 1024 5.30× 1024 1.66× 1025 5.18× 1025

0.7 1.68× 1024 5.30× 1024 1.66× 1025 5.18× 1025 9.64× 1025 2.90× 1026 8.73× 1026 2.61× 1027

0.8 9.64× 1025 2.90× 1026 8.73× 1026 2.61× 1027 3.41× 1027 9.91× 1027 2.87× 1028 8.28× 1028

0.9 3.41× 1027 9.91× 1027 2.87× 1028 8.28× 1028 8.25× 1028 2.32× 1029 6.51× 1029 1.82× 1030

1.0 8.25× 1028 2.32× 1029 6.51× 1029 1.82× 1030 8.25× 1028 2.32× 1029 6.51× 1029 1.82× 1030

Figure 1 : Approximate solution curves using classical FRPSM, example 4.1

G(t) =

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

1G(t) for t ∈ [0, 0.1],
2G(t) for t ∈ [0.1, 0.2],
3G(t) for t ∈ [0.2, 0.3],
4G(t) for t ∈ [0.3, 0.4],
5G(t) for t ∈ [0.4, 0.5],
6G(t) for t ∈ [0.5, 0.6],
7G(t) for t ∈ [0.6, 0.7],
8G(t) for t ∈ [0.7, 0.8],
9G(t) for t ∈ [0.8, 0.9],
10G(t) for t ∈ [0.9, 1],

W (t) =

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

1W (t) for t ∈ [0, 0.1],
2W (t) for t ∈ [0.1, 0.2],
3W (t) for t ∈ [0.2, 0.3],
4W (t) for t ∈ [0.3, 0.4],
5W (t) for t ∈ [0.4, 0.5],
6W (t) for t ∈ [0.5, 0.6],
7W (t) for t ∈ [0.6, 0.7],
8W (t) for t ∈ [0.7, 0.8],
9W (t) for t ∈ [0.8, 0.9],
10W (t) for t ∈ [0.9, 1],

and the Nth truncated FPS solution has the form:



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2947

G(t) =

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

1GN (t) for t ∈ [0, 0.1],
2GN (t) for t ∈ [0.1, 0.2],
3GN (t) for t ∈ [0.2, 0.3],
4GN (t) for t ∈ [0.3, 0.4],
5GN (t) for t ∈ [0.4, 0.5],
6GN (t) for t ∈ [0.5, 0.6],
7GN (t) for t ∈ [0.6, 0.7],
8GN (t) for t ∈ [0.7, 0.8],
9GN (t) for t ∈ [0.8, 0.9],
10GN (t) for t ∈ [0.9, 1],

W (t) =

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

1WN (t) for t ∈ [0, 0.1],
2WN (t) for t ∈ [0.1, 0.2],
3WN (t) for t ∈ [0.2, 0.3],
4WN (t) for t ∈ [0.3, 0.4],
5WN (t) for t ∈ [0.4, 0.5],
6WN (t) for t ∈ [0.5, 0.6],
7WN (t) for t ∈ [0.6, 0.7],
8W (t) for t ∈ [0.7, 0.8],
9WN (t) for t ∈ [0.8, 0.9],
10WN (t) for t ∈ [0.9, 1],

nGN (t) =
∑N

v=0
avn(t−tn−1)vβ

Γ(vβ+1) and nWN (t) =
∑N

v=0
bvn(t−tn−1)vβ

Γ(vβ+1) .

For N = 3, the results are as follows:

1G3(t) = 1− 9404t2β

Γ(2β + 1)
+

903544t3β

Γ(3β + 1)
+

94tβ

Γ(β + 1)
,

2G3(t) = 1 +
94(0.1)β

Γ(β + 1)
− 9404(0.1)2β

Γ(2β + 1)
+

903544(0.1)3β

Γ(3β + 1)
− 3.277291397699514(t− 0.1)β

Γ(β + 1)

+
3.5022880095564233 (t− 0.1)2β

Γ(2β + 1)
+

286.01572342177724(t− 0.1)3β

Γ(3β + 1)
,

3G3(t) =1 +
90.72270860230049(0.1)β

Γ(β + 1)
− 9400.497711990443(0.1)2β

Γ(2β + 1)

+
903830.0157234218(0.1)3β

Γ(3β + 1)
− 2.709793687108408(t− 0.2)β

Γ(β + 1)

+
5.418576468893888(t− 0.2)2β

Γ(2β + 1)
− 10.74010602678672(t− 0.2)3β

Γ(3β + 1)
,



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2948

4G3(t) =1− 2.709793687108408(0.1)β

Γ(β + 1)
+

90.72270860230049(0.1)β

Γ(β + 1)

+
5.418576468893888(0.1)2β

Γ(2β + 1)
− 9400.497711990443(0.1)2β

Γ(2β + 1)

− 10.74010602678672(0.1)3β

Γ(3β + 1)
+

903830.0157234218 (0.1)3β

Γ(3β + 1)

− 2.218600227456377(t− 0.3)β

Γ(β + 1)
+

4.437200120106484(t− 0.3)2β

Γ(2β + 1)

− 8.874368098810919(t− 0.3)3β

Γ(3β + 1)
,

5G3(t) =1− 2.709793687108408(0.1)β

Γ(β + 1)
+

90.72270860230049(0.1)β

Γ(β + 1)

− 2.218600227456377(0.1)β

Γ(β + 1)
+

5.418576468893888(0.1)2β

Γ(2β + 1)

− 9400.497711990443(0.1)2β

Γ(2β + 1)
+

4.437200120106484(0.1)2β

Γ(2β + 1)

− 10.74010602678672(0.1)3β

Γ(3β + 1)
+

903830.0157234218(0.1)3β

Γ(3β + 1)

− 8.874368098810919(0.1)3β

Γ(3β + 1)
− 1.816436237919282(t− 0.4)β

Γ(β + 1)

+
3.632872475726769(t− 0.4)2β

Γ(2β + 1)
− 7.265744940721219(t− 0.4)3β

Γ(3β + 1)
,

6G3(t) =1− 4.52622992502769(0.1)β

Γ(β + 1)
+

90.72270860230049(0.1)β

Γ(β + 1)

− 2.218600227456377(0.1)β

Γ(β + 1)
+

9.051448944620656(0.1)2β

Γ(2β + 1)

− 9400.497711990443(0.1)2β

Γ(2β + 1)
+

4.437200120106484(0.1)2β

Γ(2β + 1)

− 18.00585096750794(0.1)3β

Γ(3β + 1)
+

903830.0157234218(0.1)3β

Γ(3β + 1)

− 8.874368098810919(0.1)3β

Γ(3β + 1)
− 1.4871722089907775(t− 0.5)β

Γ(β + 1)

+
2.9743444179828282(t− 0.5)2β

Γ(2β + 1)
− 5.948688836087882(t− 0.5)3β

Γ(3β + 1)
,



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2949

7G3(t) =1− 6.013402134018468(0.1)β

Γ(β + 1)
+

90.72270860230049(0.1)β

Γ(β + 1)

− 2.218600227456377(0.1)β

Γ(β + 1)
+

12.025793362603483(0.1)2.0β

Γ(2.0β + 1)

− 9400.497711990443(0.1)2.0β

Γ(2.0β + 1)
+

4.437200120106484(0.1)2.0β

Γ(2.0β + 1)

− 23.95453980359582(0.1)3.0β

Γ(3.0β + 1.0)
+

903830.0157234218(0.1)3.0β

Γ(3.0β + 1)

− 8.874368098810919(0.1)3.0β

Γ(3.0β + 1)
− 1.2175936226236566(t− 0.6)β

Γ(β + 1.0)

+
2.4351872452475702(t− 0.6)2.0β

Γ(2.0β + 1)
− 4.870374490519822 (t− 0.6)3.0β

Γ(3.0β + 1)
,

8G3(t) =1− 7.230995756642124(0.1)β

Γ(β + 1)
+

90.72270860230049(0.1)β

Γ(β + 1)

− 2.218600227456377(0.1)β

Γ(β + 1)
+

14.460980607851054(0.1)2β

Γ(2β + 1)

− 9400.497711990443(0.1)2β

Γ(2β + 1)
+

4.437200120106484(0.1)2β

Γ(2β + 1)

− 28.824914294115644(0.1)3β

Γ(3β + 1)
+

903830.0157234218 · (0.1)3β

Γ(3β + 1)

− 8.874368098810919(0.1)3β

Γ(3β + 1)
− 0.9968813435936204(t− 0.7)β

Γ(β + 1)

+
1.9937626871880143(t− 0.7)2β

Γ(2β + 1)
− 3.9875253744502994(t− 0.7)3β

Γ(3β + 1)
,

9G3(t) =1− 7.230995756642124 (0.1)β

Γ(β + 1)
+

90.72270860230049(0.1)β

Γ(β + 1)
− 2.218600227456377(0.1)β

Γ(β + 1)

− 0.9968813435936204(0.1)β

Γ(β + 1)
+

14.460980607851054(0.1)2β

Γ(2β + 1)
− 9400.497711990443(0.1)2β

Γ(2β + 1)

+
4.437200120106484(0.1)2β

Γ(2β + 1)
+

1.9937626871880143(0.1)2β

Γ(2β + 1)
− 28.824914294115644(0.1)3β

Γ(3β + 1)

+
903830.0157234218(0.1)3β

Γ(3β + 1)
− 8.874368098810919(0.1)3β

Γ(3β + 1)
− 3.9875253744502994(0.1)3β

Γ(3β + 1)

− 0.8161774131697912(t− 0.8)β

Γ(β + 1)
+

1.6323548263398586(t− 0.8)2β

Γ(2β + 1)
− 3.2647096527062445(t− 0.8)3β

Γ(3β + 1)
.



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2950

10G3(t) =1− 8.047173169811915 · (0.1)β

Γ(β + 1)
+

90.72270860230049(0.1)β

Γ(β + 1)
− 2.218600227456377(0.1)β

Γ(β + 1)

− 0.9968813435936204(0.1)β

Γ(β + 1)
+

16.093335434190912(0.1)2β

Γ(2β + 1)
− 9400.497711990443(0.1)2β

Γ(2β + 1)

+
4.437200120106484(0.1)2β

Γ(2β + 1)
+

1.9937626871880143(0.1)2β

Γ(2β + 1)
− 32.08962394682189(0.1)3β

Γ(3β + 1)

+
903830.0157234218(0.1)3β

Γ(3β + 1)
− 8.874368098810919(0.1)3β

Γ(3β + 1)
− 3.9875253744502994(0.1)3β

Γ(3β + 1)

− 0.66822954813(t− 0.9)β

Γ(β + 1)
+

1.33645909626(t− 0.9)2β

Γ(2β + 1)
− 2.67291819253(t− 0.9)3β

Γ(3β + 1)
.

1W3(t) = 1 +
9412 t2β

Γ(2β + 1)
− 903560 t3β

Γ(3β + 1)
− 98 tβ

Γ(β + 1)
,

2W3(t) =1− 98(0.1)β

Γ(β + 1)
+

9412(0.1)2β

Γ(2β + 1)

− 903560(0.1)3β

Γ(3β + 1)
+

0.0023683853879674643(t− 0.1)β

Γ(β + 1)

+
3.04755801506667(t− 0.1)2β

Γ(2β + 1)
− 299.1154154710234(t− 0.1)3β

Γ(3β + 1)
,

3W3(t) =1− 97.99763161461203(0.1)β

Γ(β + 1)
+

9415.047558015067(0.1)2β

Γ(2β + 1)

− 903859.115415471(0.1)3β

Γ(3β + 1)
+

0.02851350296616295(t− 0.2)β

Γ(β + 1)

− 0.056016100609398656(t− 0.2)2β

Γ(2β + 1)
+

0.0149852902177372(t− 0.2)3β

Γ(3β + 1)
,

4W3(t) =1 +
0.02851350296616295(0.1)β

Γ(β + 1)
− 97.99763161461203(0.1)β

Γ(β + 1)

− 0.056016100609398656(0.1)2β

Γ(2β + 1)
+

9415.047558015067(0.1)2β

Γ(2β + 1)

+
0.0149852902177372(0.1)3β

Γ(3β + 1)
− 903859.115415471(0.1)3β

Γ(3β + 1)

+
0.023353683080527432(t− 0.3)β

Γ(β + 1)
− 0.04670703135478371(t− 0.3)2β

Γ(2β + 1)

+
0.0933819213075534(t− 0.3)3β

Γ(3β + 1)
,



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2951

5W3(t) =1 +
0.02851350296616295(0.1)β

Γ(β + 1)
− 97.99763161461203(0.1)β

Γ(β + 1)

+
0.023353683080527432(0.1)β

Γ(β + 1)
− 0.056016100609398656(0.1)2β

Γ(2β + 1)

+
9415.047558015067(0.1)2β

Γ(2β + 1)
− 0.04670703135478371(0.1)2β

Γ(2β + 1)

+
0.0149852902177372(0.1)3β

Γ(3β + 1)
− 903859.115415471(0.1)3β

Γ(3β + 1)

+
0.0933819213075534(0.1)3β

Γ(3β + 1)
+

0.01912038145060513(t− 0.4)β

Γ(β + 1)

− 0.03824076278941546(t− 0.4)2β

Γ(2β + 1)
+

0.07648151484653454(t− 0.4)3β

Γ(3β + 1)
,

6W3(t) =1 +
0.04763388441676808(0.1)β

Γ(β + 1)
− 97.99763161461203(0.1)β

Γ(β + 1)
+

0.023353683080527432(0.1)β

Γ(β + 1)

− 0.09425686339881412(0.1)2β

Γ(2β + 1)
+

9415.047558015067(0.1)2β

Γ(2β + 1)
− 0.04670703135478371(0.1)2β

Γ(2β + 1)

+
0.09146680506427174(0.1)3β

Γ(3β + 1)
− 903859.115415471(0.1)3β

Γ(3β + 1)
+

0.0933819213075534(0.1)3β

Γ(3β + 1)

+
0.015654444305179482(t− 0.5)β

Γ(β + 1)
− 0.03130888861163217(t− 0.5)2β

Γ(2β + 1)
+

0.06261777734550833(t− 0.5)3β

Γ(3β + 1)
,

7W3(t) =1 +
0.06328832872194756(0.1)β

Γ(β + 1)
− 97.99763161461203(0.1)β

Γ(β + 1)
+

0.023353683080527432(0.1)β

Γ(β + 1)

− 0.12556575201044629(0.1)2β

Γ(2β + 1)
+

9415.047558015067(0.1)2β

Γ(2β + 1)
− 0.04670703135478371(0.1)2β

Γ(2β + 1)

+
0.15408458240978007(0.1)3β

Γ(3β + 1)
− 903859.115415471(0.1)3β

Γ(3β + 1)
+

0.0933819213075534(0.1)3β

Γ(3β + 1)

+
0.012816774974988565(t− 0.6)β

Γ(β + 1)
− 0.025633549950234258(t− 0.6)2β

Γ(2β + 1)
+

0.05126709992515507(t− 0.6)3β

Γ(3β + 1)

8W3(t) =1 +
0.07610510369693613(0.1)β

Γ(β + 1)
− 97.99763161461203(0.1)β

Γ(β + 1)
+

0.023353683080527432(0.1)β

Γ(β + 1)

− 0.15119930196068054(0.1)2β

Γ(2β + 1)
+

9415.047558015067(0.1)2β

Γ(2β + 1)
− 0.04670703135478371(0.1)2β

Γ(2β + 1)

+
0.20535168233493514(0.1)3β

Γ(3β + 1)
− 903859.115415471(0.1)3β

Γ(3β + 1)
+

0.0933819213075534(0.1)3β

Γ(3β + 1)

+
0.010493487827309411(t− 0.7)β

Γ(β + 1)
− 0.020986975655392648(t− 0.7)2β

Γ(2β + 1)
+

0.041973951385060104(t− 0.7)3β

Γ(3β + 1)
,



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2952

Figure 2 :Approximate solution curves using MS-FRPSM, example 4.1

9W3(t) =1 +
0.07610510369693613(0.1)β

Γ(β + 1)
− 97.99763161461203(0.1)β

Γ(β + 1)
+

0.023353683080527432(0.1)β

Γ(β + 1)
+

0.010493487827309411(0.1)β

Γ(β + 1)

− 0.15119930196068054(0.1)2β

Γ(2β + 1)
+

9415.047558015067(0.1)2β

Γ(2β + 1)
− 0.04670703135478371(0.1)2β

Γ(2β + 1)
− 0.020986975655392648(0.1)2β

Γ(2β + 1)

+
0.20535168233493514(0.1)3β

Γ(3β + 1)
− 903859.115415471(0.1)3β

Γ(3β + 1)
+

0.0933819213075534(0.1)3β

Γ(3β + 1)
+

0.041973951385060104(0.1)3β

Γ(3β + 1)

+
0.008591341191263868(t− 0.8)β

Γ(β + 1)
− 0.01718268238280407(t− 0.8)2β

Γ(2β + 1)
+

0.03436536479212293(t− 0.8)3β

Γ(3β + 1)
+ 1,

10W3(t) =1 +
0.0846964448882(0.1)β

Γ(β + 1)
− 97.99763161461203(0.1)β

Γ(β + 1)
+

0.023353683080527432(0.1)β

Γ(β + 1)
+

0.010493487827309411(0.1)β

Γ(β + 1)

− 0.16838198434348461(0.1)2β

Γ(2β + 1)
+

9415.047558015067(0.1)2β

Γ(2β + 1)
− 0.04670703135478371(0.1)2β

Γ(2β + 1)
− 0.020986975655392648(0.1)2β

Γ(2β + 1)

+
0.23971704712705807(0.1)3β

Γ(3β + 1)
− 903859.115415471(0.1)3β

Γ(3β + 1)
+

0.0933819213075534(0.1)3β

Γ(3β + 1)
+

0.041973951385060104(0.1)3β

Γ(3β + 1)

+
0.00703399524347198(t− 0.9)β

Γ(β + 1)
− 0.014067990487040993(t− 0.9)2β

Γ(2β + 1)
+

0.02813598098340719(t− 0.9)3β

Γ(3β + 1)
.

The MS-FRPS is carried out for N = 30 and some values of the approximate numerical
solutions, together with the absolute errors for β = 1, are presented in Table 3. The
approximate solution curves and values of the residual errors for different values of β are
given in Figure 2 and Table 4, respectively.

Obviously, the MS-FRPSM gives more accurate results than those found with classical
FRPSM in Tables 1 and 2, and Figure 1. However, the residual error is still large, which
can be decreased by performing more iterations.

Example 4.2: Consider the following nonhomogeneous fractional stiff system:

C
t D

β
0G(t) = −G(t)− 15W (t)− 15e−t, t ∈ [0, 1], 0 < β ≤ 1,

C
t D

β
0W (t) = 15G(t)−W (t)− 15e−t,

G(0) = 1, W (0) = 1.

(12)

The exact solution for the system (12) when β = 1 is:

G(t) = e−t, W (t) = e−t.



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2953

Table 3 Numerical results for MS-FRPSM for example 4.1

ti G(t) G30(t), β = 1 |G(t)−G30(t)|, β = 1 G30(t), β = 0.99 G30(t), β = 0.98 G30(t), β = 0.97

0.0 1.00000 1.00000 0.00000 1.00000 1.00000 1.00000
0.1 1.65481 1.65454 0.000269072 1.64314 1.62541 1.57125
0.2 1.35490 1.35490 1.0734 ×10−7 1.33634 1.31133 1.24976
0.3 1.10930 1.10930 3.65752 ×10−11 1.08487 1.05391 0.986293
0.4 0.90822 0.90822 1.72085×10−13 0.87899 0.84315 0.77058
0.5 0.74359 0.74359 1.94289 ×10−13 0.71043 0.67059 0.59397
0.6 0.60880 0.60880 5.3435 ×10−14 0.57242 0.52931 0.44938
0.7 0.49844 0.49844 3.57769 ×10−13 0.45943 0.41365 0.33099
0.8 0.40809 0.40809 1.32228×10−13 0.36692 0.31895 0.23407
0.9 0.33411 0.33411 2.32314 ×10−14 0.29118 0.24141 0.15471
1.0 0.27355 0.27355 2.37588 ×10−14 0.22917 0.17793 0.08974

ti W (t) W30(t), β = 1 |W (t)−W30(t)|, β = 1 W30(t), β = 0.99 W30(t), β = 0.98 W30(t), β = 0.97

0.0 1.00000 1.00000 0.00000 1.00000 1.00000 1.00000
0.1 -0.0173506 -0.0170816 0.000269072 -0.0143377 -0.00539508 0.0398389
0.2 -0.0142621 -0.0142620 1.07339× 10−7 -0.0114419 -0.00241976 0.0429069
0.3 -0.0116768 -0.0116768 3.64694× 10−11 -0.00879496 0.00028984 0.0456802
0.4 -0.00956019 -0.00956019 5.91012× 10−13 -0.00662777 0.00250836 0.0479508
0.5 -0.00782722 -0.00782722 3.56415× 10−13 -0.00485343 0.00432473 0.0498099
0.6 -0.00640839 -0.00640839 6.70992× 10−13 -0.00340072 0.00581185 0.0513319
0.7 -0.00524674 -0.00524674 5.19546× 10−13 -0.00221135 0.0070294 0.0525781
0.8 -0.00429567 -0.00429567 2.34957× 10−13 -0.00123757 0.00802624 0.0535984
0.9 -0.00351700 -0.00351700 1.33380× 10−13 -0.000440303 0.00884239 0.0544337
1.0 -0.00287947 -0.00287947 1.06159× 10−13 0.000212441 0.0095106 0.0551176

Table 4 Some numerical values of the residual function for MS-FRPSM, Example 4.1

ti |ResidG,30(t)| |ResidW,30(t)|
β = 1 β = 0.99 β = 0.98 β = 1 β = 0.99 β = 0.98

0.1 3.27729 3.00522 2.13794 0.00236839 0.25239 1.10209
0.2 2.70979 2.42332 1.5412 0.0285135 0.22648 1.07661
0.3 2.2186 1.9204 1.02637 0.0233537 0.231763 1.08202
0.4 1.81644 1.50863 0.604851 0.0191204 0.236098 1.08646
0.5 1.48717 1.17151 0.259741 0.0156544 0.239646 1.09009
0.6 1.21759 0.895491 0.0228115 0.0128168 0.242552 1.09306
0.7 0.996881 0.669509 0.254146 0.0104935 0.24493 1.0955
0.8 0.816177 0.484491 0.443547 0.00859134 0.246878 1.09749
0.9 0.66823 0.333011 0.598615 0.007034 0.248472 1.09912
1.0 0.5471 0.208989 0.725574 0.00575895 0.249778 1.10046

To apply the MS-FRPSM, we divide the interval [0,1] into 10 subintervals, so we have
step size h= 0.1. On each subinterval, we carry out 10 iterations of the FRPSM. Accurate
results for stiff system in (12) are clear from the absolute errors listed in Table 5 and the
residual errors in Table 6. Moreover, the approximate solution curves are displayed in
Figure 3, while the behavior of the 10th residual functions is shown in Figure 4.



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2954

Table 5 Numerical results for MS-FRPSM for example 4.2

ti G(t) G10(t), β = 1 |G(t)−G10(t)|, β = 1 G10(t), β = 0.95 G10(t), β = 0.9 G10(t), β = 0.85

0.0 1 1 0 1 1 1
0.1 0.904837 0.904837 1.11022 ×10−16 0.892109 0.878096 0.862774
0.2 0.818731 0.818731 2.22045 ×10−16 0.794485 0.767793 0.738607
0.3 0.740818 0.740818 2.22045 ×10−16 0.706152 0.667986 0.626256
0.4 0.67032 0.67032 2.22045 ×10−16 0.626224 0.577678 0.524597
0.5 0.606531 0.606531 1.11022 ×10−16 0.553903 0.495963 0.432612
0.6 0.548812 0.548812 1.11022 ×10−16 0.488463 0.422025 0.349380
0.7 0.496585 0.496585 5.55112 ×10−17 0.429252 0.355123 0.274069
0.8 0.449329 0.449329 0 0.375675 0.294587 0.205925
0.9 0.40657 0.40657 1.11022 ×10−16 0.327196 0.239812 0.144265
1.0 0.367879 0.367879 0 0.283331 0.19025 0.0884732

ti W (t) W30(t), β = 1 |W (t)−W10(t)|, β = 1 W10(t), β = 0.95 W10(t), β = 0.9 W10(t), β = 0.85

0.0 1.000000 1.000000 0.000000 1.000000 1.000000 1.000000
0.1 0.904837 0.904837 1.11022× 10−16 0.892109 0.878096 0.862774
0.2 0.818731 0.818731 1.11022× 10−16 0.794485 0.767793 0.738607
0.3 0.740818 0.740818 1.11022× 10−16 0.706152 0.667986 0.626256
0.4 0.670320 0.670320 1.11022× 10−16 0.626224 0.577678 0.524597
0.5 0.606531 0.606531 1.11022× 10−16 0.553903 0.495963 0.432612
0.6 0.548812 0.548812 0.000000 0.488463 0.422025 0.349380
0.7 0.496585 0.496585 1.66533× 10−16 0.429252 0.355123 0.274069
0.8 0.449329 0.449329 1.66533× 10−16 0.375675 0.294587 0.205925
0.9 0.406570 0.406570 5.55112× 10−17 0.327196 0.239812 0.144265
1.0 0.367879 0.367879 0.000000 0.283331 0.190250 0.0884732

Table 6 Some numerical values of the residual function for MS-FRPSM, example 4.2

ti |ResidG,10(t)| |ResidW,10(t)|
β = 0.95 β = 0.90 β = 0.85 β = 0.95 β = 0.90 β = 0.85

0.1 0.00018 0.00018 0.00017 0.00011 0.00013 0.00015
0.2 0.00015 0.00015 0.00014 0.00012 0.00015 0.00020
0.3 0.00012 0.00011 0.00010 0.00012 0.00018 0.00023
0.4 0.00010 0.00009 0.00008 0.00013 0.00020 0.00027
0.5 0.00008 0.00007 0.00006 0.00013 0.00022 0.00030
0.6 0.00006 0.00005 0.00004 0.00014 0.00023 0.00033
0.7 0.00005 0.00004 0.00003 0.00014 0.00025 0.00036
0.8 0.00004 0.00003 0.00002 0.00015 0.00026 0.00039
0.9 0.00003 0.00002 0.00001 0.00016 0.00027 0.00041
1.0 0.00002 0.00001 7.2 ×10−6 0.00016 0.00029 0.00043



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2955

Figure 3 : Approximate solution curves using MS-FRPSM, example 4.2

Figure 4 : The 10th residual function using MS-FRPSM, example 4.2

Example 4.3: Consider the following fractional stiff system:

C
t D

β
0G(t) = −2G(t) +W (t) + 2 sin(t), t ∈ [0, 5], 0 < β ≤ 1, (13)

C
t D

β
0W (t) = −3G(t) + 2 (W (t) + sin(t)− cos(t)),

with initial conditions
G(0) = 2, W (0) = 3. (14)

The exact solution of (13) when β = 1 is G(t) = 2e−t + sin t, W (t) = 2e−t + cos t. We
solve the fractional stiff system in (13) using classical FRPSM and the MS-FRPSM with
step size h = 1. In both methods, we compute the 10th truncated FPS solution, which is
denoted by G10,10(t) and W10,10(t) for MS-FRPSM, and by G10(t) and W10(t) for classical
FRPSM. A comparison between these two methods is carried out as follows: Table 7
compares the absolute errors for β = 1, while Table 8 and Table 9 represent comparisons
between the residual errors for β = 0.9, β = 0.8, and β = 0.7. It is noticeable from these
tables that multistep schemes contribute to reducing the error of approximations even
though it is still high. The graphs of the 10th residual functions are shown in Figure 5
and Figure 6. On the other hand, a comparison between the curves of the 10th FPS using
both methods for β = 0.9, β = 0.8, and β = 0.7 are presented in Figure 7 and Figure 8.



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2956

Table 7 A comparison between absolute errors using FRPSM and MS-FRPSM when
β = 1, Example 4.3

MS-FRPSM FRPSM MS-FRPSM FRPSM

ti |G(t)−G10,10(t)| |G(t)−G10(t)| |W (t)−W10,10(t)| |W (t)−W10(t)|
0.0 0.000000 0.000000 0.000000 0.000000
0.5 3.56963× 10−11 3.56963× 10−11 2.29745× 10−11 2.29745× 10−11

1.0 7.11208× 10−8 7.11208× 10−8 4.41523× 10−8 4.41523× 10−8

1.5 2.90986× 10−8 5.98460× 10−6 1.53817× 10−8 3.58103× 10−6

2.0 2.31792× 10−8 0.000137827 8.38990× 10−8 0.000079445
2.5 4.17414× 10−8 0.00156031 2.18289× 10−7 0.00086600
3.0 1.23302× 10−7 0.01127010 4.23894× 10−7 0.00602104
3.5 2.31435× 10−7 0.05969050 7.27017× 10−7 0.03069310
4.0 4.21251× 10−7 0.25188200 1.21487× 10−6 0.12466700
4.5 6.69059× 10−7 0.89344500 1.97749× 10−6 0.42573800
5.0 1.10148× 10−6 2.76316000 3.22386× 10−6 1.26819000

Table 8 A comparison between residual errors |ResidG| using MS-FRPSM and FRPSM,
Example 4.3

ti MS-FRPSM FRPSM

β = 0.9 β = 0.8 β = 0.7 β = 0.9 β = 0.8 β = 0.7

0.5 0.413857 0.494220 0.566087 0.413857 0.494224 0.566087
1.0 0.231610 0.226850 0.181151 0.231610 0.226850 0.181151
1.5 0.425311 0.412688 0.331603 0.246721 0.084573 0.080431
2.0 0.608996 0.469273 0.283837 0.388434 0.132785 0.034489
2.5 0.788556 0.530022 0.210149 0.626232 0.442167 0.472802
3.0 0.864512 0.599768 0.327769 0.953735 1.09011 1.65211
3.5 0.638268 0.241984 0.148857 1.441670 2.29235 4.00186
4.0 0.424538 0.165570 0.053341 2.406320 4.70177 8.68621
4.5 0.082856 0.369207 0.585395 4.77180 9.92627 18.0669
5.0 0.431016 0.556547 0.628991 10.74120 21.3298 36.3969



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2957

Table 9 A comparison between residual errors |ResidW | using FRPSM and MS-FRPSM,
Example 4.3

ti MS-FRPSM FRPSM

β = 0.9 β = 0.8 β = 0.7 β = 0.9 β = 0.8 β = 0.7

0.5 7.61053 7.43003 7.23856 7.61053 7.43003 7.23856
1.0 4.35363 4.25856 4.15649 4.35363 4.25856 4.15649
1.5 2.27806 2.08497 1.83811 2.28767 2.16498 2.08242
2.0 1.68566 1.46167 1.19799 1.51504 1.26703 1.15864
2.5 1.85531 1.47128 0.98933 1.56928 1.26046 1.27452
3.0 1.86966 1.45236 1.00075 1.78177 1.79020 2.43064
3.5 0.92400 0.297301 0.34279 1.93607 3.02567 5.35909
4.0 0.22727 0.595633 0.89744 2.86553 6.44237 12.6162
4.5 1.17682 1.429970 1.54596 7.21321 16.2054 30.3349
5.0 0.812702 0.773654 0.63014 20.9921 41.3743 70.2115

Figure 5 : Curves of the 10th residual functions for G(t) of system (13) using
MS-FRPSM (right) and FRPSM (left)

Figure 6 :Curves of the 10th residual functions for W(t) of system (13) using
MS-FRPSM (right) and FRPSM (left)



M. Al Zurayqat, S. Hasan / Eur. J. Pure Appl. Math, 17 (4) (2024), 2939-2961 2958

Figure 7 : The 10th FPS of G(t) for system (13) using MS-FRPSM (right) and FRPSM
(left)

Figure 8 : The 10th FPS of W(t) for system (13) using MS-FRPSM (right) and FRPSM
(left)

5. Conclusions

In this paper, we introduced a modified algorithm, namely, the MS-FRPSM, in an
attempt to obtain approximate solutions of fractional stiff systems. Numerical examples
were carried out to assess the applicability of the improved technique for solving these
systems. We compared the exact solution for integer order stiff systems to the N -th
truncated FPS, in addition to comparing the residual errors for fractional orders obtained
using classical FRPSM and MS-FRPSM.
The comparison reveals that MS-FRPSM reduced both absolute and residual errors. More
iterations and a smaller step size lead to higher accuracy. Moreover, MS-FRPSM increased
the length of the convergence interval. Additionally, it does not require large computer
memory and provides accurate numerical results in less time.
For future works, it is recommended to apply the MS-FRPSM for more types of equations
with different fractional operators. It may be combined with some decomposition methods
such as Adomian decomposition method to solve nonlinear FDEs.



REFERENCES 2959

References

[1] M Adel, MM Khader, Hijaz Ahmad, and TA Assiri. Approximate analytical solutions
for the blood ethanol concentration system and predator-prey equations by using
variational iteration method. Aims Math, 8(8):19083–19096, 2023.

[2] Mohamed Adel, Mohamed M Khader, Taghreed A Assiri, and Wajdi Kallel. Numeri-
cal simulation for covid-19 model using a multidomain spectral relaxation technique.
Symmetry, 15(4):931, 2023.

[3] OA Akinfenwa, B Akinnukawe, and SB Mudasiru. A family of continuous third
derivative block methods for solving stiff systems of first order ordinary differential
equations. Journal of the Nigerian Mathematical Society, 34(2):160–168, 2015.

[4] Hossein Aminikhah and Milad Hemmatnezhad. An effective modification of the homo-
topy perturbation method for stiff systems of ordinary differential equations. Applied
Mathematics Letters, 24(9):1502–1508, 2011.

[5] Mehmet Tarik Atay and Okan Kilic. The semianalytical solutions for stiff systems
of ordinary differential equations by using variational iteration method and modi-
fied variational iteration method with comparison to exact solutions. Mathematical
Problems in Engineering, 2013(1):143915, 2013.

[6] Shalashilin Vladimir Ivanovich and Kuznetsov Evgenĭı Borisovich. Parametric contin-
uation and optimal parametrization in applied mathematics and mechanics. Springer
Science & Business Media, 2013.

[7] S. Cifani and E. R. Jakobsen. Entropy solution theory for fractional degenerate
convection–diffusion equations. Annales de l’Institut Henri Poincaré C, 28(3):413–
441, 2011.

[8] Charles Francis Curtiss and Joseph O Hirschfelder. Integration of stiff equations.
Proceedings of the national academy of sciences, 38(3):235–243, 1952.

[9] Ahmad El-Ajou, Omar Abu Arqub, and Mohammed Al-Smadi. A general form of
the generalized taylor’s formula with some applications. Applied Mathematics and
Computation, 256:851–859, 2015.

[10] A. Freihet, S. Hasan, M. Al-Smadi, M. Gaith, and S. Momani. Construction of
fractional power series solutions to fractional stiff system using residual functions
algorithm. Advances in Difference Equations, 2019(1):1–15, 2019.

[11] S. Hasan, M. Al-Smadi, S. Momani, and O. A. Arqub. Residual power series approach
for solving linear fractional swift-hohenberg problems. In Mathematical Methods and
Modelling in Applied Sciences, pages 33–43. Springer International Publishing, 2020.



REFERENCES 2960

[12] Shatha Hasan, Mohammed Al-Smadi, Hemen Dutta, Shaher Momani, and Samir
Hadid. Multi-step reproducing kernel algorithm for solving caputo–fabrizio fractional
stiff models arising in electric circuits. Soft Computing, 26(8):3713–3727, 2022.

[13] G. M. Ismail, H. R. Abdl-Rahim, H. Ahmad, and Y. M. Chu. Fractional residual
power series method for the analytical and approximate studies of fractional physical
phenomena. Open Physics, 18(1):799–805, 2020.

[14] Ali Jameel, NR Anakira, AK Alomari, Noraziah H Man, et al. Solution and analysis
of the fuzzy volterra integral equations via homotopy analysis method. Computer
Modeling in Engineering & Sciences, 127(3):875–899, 2021.

[15] Ali Fareed Jameel, Nidal Ratib Anakira, AK Alomari, DM Alsharo, and Azizan Saa-
ban. New semi-analytical method for solving two point nth order fuzzy boundary
value problem. International Journal of Mathematical Modelling and Numerical Op-
timisation, 9(1):12–31, 2019.

[16] M. Khader and M. H. DarAssi. Residual power series method for solving non-
linear reaction-diffusion-convection problems. Boletim da Sociedade Paranaense de
Matemática, 39(3):177–188, 2021.

[17] Z. Korpinar, M. Inc, E. Hınçal, and D. Baleanu. Residual power series algorithm for
fractional cancer tumor models. Alexandria Engineering Journal, 59(3):1405–1412,
2020.

[18] J.A.T. Machado. Entropy analysis of integer and fractional dynamical systems. Non-
linear Dynamics, 62:371–378, 2010.

[19] I. Podlubny. Fractional Differential Equations. Academic Press, San Diego, 1999.

[20] M. Qayyum and Q. Fatima. Solutions of stiff systems of ordinary differential equations
using residual power series method. Journal of Mathematics, 2022.

[21] T. R. Rao. Application of residual power series method to time fractional gas dynam-
ics equations. In Journal of Physics: Conference Series, volume 1139, page 012007.
IOP Publishing, 2018.

[22] DG Yakubu and S Markus. The efficiency of second derivative multistep methods
for the numerical integration of stiff systems. Journal of the Nigerian mathematical
Society, 35(1):107–127, 2016.

[23] Feras Yousef, Osama Alkam, and Ines Saker. The dynamics of new motion styles in
the time-dependent four-body problem: weaving periodic solutions. The European
Physical Journal Plus, 135:1–10, 2020.

[24] Feras Yousef, Billel Semmar, and Kamal Al Nasr. Dynamics and simulations of dis-
cretized caputo-conformable fractional-order lotka–volterra models. Nonlinear Engi-
neering, 11(1):100–111, 2022.



REFERENCES 2961

[25] Feras Yousef, Billel Semmar, and Kamal Al Nasr. Incommensurate conformable-
type three-dimensional lotka–volterra model: discretization, stability, and bifurcation.
Arab Journal of Basic and Applied Sciences, 29(1):113–120, 2022.