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

Exact Solutions oF Nonlinear Delay Voltera
Integro-Differential Equations Using Modified

Homotopy Perturbation Method

Nidal Anakira1,2,∗, Ala Amourah1, Adel Almalki3, Abdullah Alsoboh4,∗,
Tala Sasa5

1 Faculty of Education and Arts, Sohar University, Sohar 3111, Oman
2 Jadara University Research Center, Jadara University, Jordan
3 Al-Gunfudah University College, Umm Al- Qura University, Mecca 21955, Saudi Arabia
4 College of Applied and Health Sciences, Al Sharqiyah University, Post Box No. 42,
Post Code No. 400, Ibra, Sultanate of Oman
5 Department of Mathematics, Faculty of Science, Private Applied Science University, Am-
man, Jordan

Abstract. This paper presents a robust and efficient approach for solving delay Volterra integro-
differential equations (DVIDEs), which model systems with memory and delay effects commonly
encountered in fields such as biology, control theory, and epidemiology. Due to their complexity,
these equations require accurate and efficient solution methods. The proposed method modifies the
traditional homotopy perturbation method (HPM) to form a new version (MHPM), integrating
it with the Laplace transformation and Padé approximants. The incorporation of the Laplace
transformation simplifies the problem by converting integro-differential equations into algebraic
equations, streamlining the solution process. To further enhance the accuracy and convergence
of the solution series, Padé approximants are employed, enabling the method to overcome the
limitations of standard perturbation techniques. This hybrid approach effectively combines the
strengths of homotopy perturbation, Laplace transformation, and Padé approximants, yielding
highly accurate solutions that closely approximate the exact ones for various nonlinear DVIDEs.
Numerical experiments and illustrative examples confirm the method’s efficiency and superior
accuracy, even for equations with complex delay terms. The results highlight the potential of this
combined approach as a powerful analytical tool for solving nonlinear delay integro-differential
equations of the Volterra type in scientific and engineering applications.

2020 Mathematics Subject Classifications: 26A33, 11Cxx, 34A08

Key Words and Phrases: Homotopy Perturbation Method (HPM), Nonlinear Delay Volterra
Integro-Differential Equations, Laplace Transform, Padé Approximants

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

Email addresses: nanakira@su.edu.om (N. Anakira), AAmourah@su.edu.om (A. Amourah),
aaamalki@uqu.edu.su (A. Almalki), abdullah.alsoboh@asu.edu.om (A. Alsoboh),
t sasa@asu.edu.jo (T. Sasa)

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



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 2 of 16

1. Introduction

Ordinary, partial, integral, fractional, fuzzy, and functional differential equations are
vital tools for modeling complex phenomena in science and engineering. Each type cap-
tures distinct dependencies and dynamic behaviors, making them essential for accurately
describing real-world systems. To solve these equations effectively, researchers have devel-
oped a wide range of analytical and numerical methods, including the homotopy pertur-
bation method (HPM), the Adomian decomposition method (ADM), and the variational
iteration method (VIM), along with various modifications. These techniques enhance the
accuracy, convergence, and applicability of solutions across diverse problems. For instance,
the optimal homotopy asymptotic Method (OHAM), originally introduced by Marinca et
al. in 2008 [1, 2], has been successfully applied to different types of differential equations,
such as Volterra integro-differential equations [3], delay differential equations [4], singular
two-point boundary value problems [5], nonlinear anharmonic oscillators [6], and fuzzy
heat equations [7]. Other widely used methods include HPM [8–11], ADM [12–15], homo-
topy analysis method (HAM) [16–19], the Collection Method [20, 21], and VIM [22–24], all
of which have become indispensable in modern applications, especially with the support
of powerful computational tools and software [25–35].

Numerical methods for solving Volterra integro-differential equations play a pivotal
role in addressing complex problems in science and engineering [31–37]. These equations,
characterized by their dependence on both the current state and past data, arise in di-
verse fields such as population dynamics, viscoelasticity, and electrical circuit analysis. A
particular subclass, DVIDEs, incorporates both integral and delay terms in addition to dif-
ferential components. These equations are especially useful for modeling systems where the
present state depends on previous states, rates of change, and delayed responses. They are
frequently used to describe processes involving hereditary effects and time lags—common
in disciplines such as engineering, biology, economics, and applied sciences. As such, the
development of efficient numerical solutions for DVIDEs is of significant interest. Among
the semi-analytical methods, the HPM and its modifications—originally introduced by Ji-
Huan He in 2003—are particularly notable for their simplicity and efficiency in solving a
wide range of linear and nonlinear differential equations. Over time, several enhancements
have been developed to improve its accuracy and applicability [38–43]. HPM combines
concepts from topology and classical perturbation theory to construct a homotopy that
smoothly transitions from a simple solvable problem to the original complex one. An
embedding parameter is introduced, and the solution is expressed as a series expansion
in terms of this parameter. Once the terms are obtained iteratively, the parameter is set
to P = 1, yielding an approximate solution. A key advantage of HPM is that it does not
require the presence of a small parameter in the original equation, making it applicable to
a wider range of problems. Additionally, HPM often converges rapidly, providing accurate
results with minimal computational effort. It has been successfully applied in fields such
as fluid dynamics, heat and mass transfer, structural analysis, and quantum mechanics.
However, the success of HPM relies on the proper construction of the homotopy, and
convergence is not guaranteed for all problems. To overcome such limitations, a modified



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 3 of 16

version of HPM (MHPM) has been proposed, integrating Laplace transforms and Padé ap-
proximants to improve solution accuracy and convergence behavior. The Modified HPM,
when combined with the Laplace transform and Padé approximants, forms a powerful
hybrid method for solving nonlinear differential equations. This approach leverages the
strengths of each component: the Modified HPM reformulates the nonlinear problem into
a more manageable form via homotopy; the Laplace transform simplifies the handling of
derivatives and initial conditions by converting the problem into the Laplace domain; and
Padé approximants enhance solution accuracy by converting truncated series into rational
functions, which often converge more rapidly and represent the solution more effectively
over a broader range. This hybrid technique is particularly effective for tackling strongly
nonlinear problems, singularities, and boundary layers—situations where standard series
solutions may diverge or converge slowly. By efficiently addressing initial and boundary
conditions and improving convergence properties, the combined use of MHPM, Laplace
transforms, and Padé approximants has shown great promise in solving problems across
various scientific and engineering domains, including fluid mechanics, heat transfer, and
nonlinear oscillations.

2. Methodology

2.1. Homotopy Perturbation Method

To explain the core idea of the Homotopy Perturbation Method (HPM), we will analyze
the following equation according to, see [8–11].

A(u)− f(r) = 0, r ∈ Ω, (1)

where A is the integral operator composed of the linear operator L and the nonlinear
operator N, f(r) is a dependent variable, and Γ is the boundary of the domain Ω. Eq. (1)
can be

L(u)−N(u)− f(r) = 0. (2)

A homotopy equation v : Ω[0, 1] → R that satisfies

H(v; p) = L(v)− L (v0) + pL (v0) + p[N(v)− f(r)] = 0, (3)

or
H(v; p) = (1− p) [L(v)− L (v0)] + p [A (v0)− F (r)] = 0, (4)

is being constructed here, r ∈ Ω, p ∈ [0, 1] is the homotopy parameter, and v0(x) is an
initial approximation of Eq. (1). It is observed that

(v; 0) = L(u)− L (v0) = 0, H(v; 1) = A(v)− F (r) = 0. (5)

The process of change p from zero to one includes the transformation of H(v; p) from
L(u) − L (v0) to A(v) − F (r) is named deformation, Furthermore, L(u) − L (v0) and



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 4 of 16

A(v) − F (r) are named homotopic. Observe that 0 ≤ p ≤ 1 is considered as a small
parameter, the solution of Eqs. (4) or (5) can be expressed as a series in p, as follows:

v = v0 + pv1 + p2v2 + p3v3 + . . . (6)

when p → 1, Eq. (2.1.4) or Eq. (5) corresponds to Eq. (3) and becomes the approximate
solution of Eq. (1). i.e.

u(x) = lim
p→1

v(x) = v1 + v2 + v3 + · · · (7)

2.2 Padè approximation

For the function u(x), the Padé approximation of order
[
L
M

]
, for more details, see

[44, 45] , can be [
L

M

]
=

PL(x)

QM (x)
,

The Padé approximant of order
[
L
M

]
to u(x) is the rational function

u[L/M ](x) =
PL(x)

QM (x)
=

b0 + b1x+ · · ·+ bmxm

1 + c1x+ · · ·+ cnxn

where PL(x) and QM (x), are two polynomials of the highest degree L and M . The
power series

u(x) =

∞∑
i=1

aix
i

The coefficients of the polynomials PL(x) and QM (x), can be obtained from

u(x)− PL(x)

QM (x)
= O

(
xL+M+1

)
(8)

When the denominator and numerator functions PL(x)
QM (x) are multiplied by a constant that

is not zero, the fractional values stay the same, so that we can set the normalization
requirement as follows:

QM(0) = 1 (9)

The polynomial associated with the functions PL(x) and QM (x) was found to have no
public factors. The coefficients of the polynomial QM (x) and PL(x) are given by

PL(t) = P0 + P1t+ P2t
2 + · · ·+ PLt

L.

QM (t) = q0 + q1t+ q2t
2 + · · ·+ qM tM . (10)



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 5 of 16

The following linear systems of coefficients can be obtained by multiplying Eq. (1) by
QM (x) considering Eq. (3).

aL+1 + aLq1 + · · ·+ aL−M+1qM = 0
aL+2 + aL+1q1 + · · ·+ aL−M+2qM = 0

·
·

aL+M + aL+M−1q1 + · · ·+ aLqM = 0

 , (11)

a0 = P0
a0 + a0q1 = P1

a2 + a1q1 + a0q2 = P2

·
·

aL + aL−1q1 + · · ·+ a0qL = PL

 . (12)

These equations will be solved using Eq. (4), It is seen as a set of linear formulas for
the unidentified q ’s. When the q ’s are identified, then (5) have an explicit formula for
the unknown p‘s, this concludes the solution to the problem. If Eqs. (4) and (5) are
non-singular, then we can solve them directly and get (6), where (6) possesses, if the lower
index on a sum exceeds the upper, the sum is replaced by zero:

[
L

M

]
=

det


aL−M+1 aL−M+2 · · · aL+1

· · · · · ·
· · · ·
aL aL+1 aL+M∑L

j=M aj−MXj
∑L

j=M−1 aj−M+1X
j . . .

∑L
j=0 ajX

j



det


aL−M+1 aL−M+2 · · · aL+1

· · · ·
· · · ·
aL aL+1 . . . aL+M

XM XM−1 . . . 1


.

3. Applications of HPM

In this section, we present two illustrative examples of nonlinear DVIDEs. These ex-
amples aim to demonstrate the effectiveness and reliability of the modified HPM procedure.

3.1. Example 1 Given the following nonlinear DVIDE [46],

u′(x) = u2
(x
2

)
+

∫ t

0
u2

(
S

2

)
dS − ex + 1, u(x) = ex, 0 ≤ x, (13)



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 6 of 16

with exact solution u(x) = ex.
To find approximate solutions to this problem, we rewrite the above problem into the
following linear and nonlinear operators:

L
(
u(x)) = u′(x),

N(u(x)) = −u2
(
S

2

)
−
∫ x

0
u2

(s
2

)
ds+

5∑
n=0

xn

n!
− 1. (14)

Then we construct the homotopy equation.

u′(x)− u0(x) + p

[
−u2

(s
2

)
−
∫ x

0
u2

(s
2

)
ds+

5∑
n=0

xn

n!
− 1

]
. (15)

Using HPM, we expand u(x) as a perturbation series in p :

v(x) = v0(x) + pv1(x) + p2v2(x) + p3v3(x) + p4v4(x) + p5v5(x) + · · · (16)

The zeroth order approximation that is obtained by setting p = 0, is as follows

u′0(x) = 0, u0(0) = 1, (17)

which has the following solutions u0(x) = 1.
Moreover, the first, second, and other orders of approximations will be obtained by

substituting the homotopy equations and matching the terms of power p to get the fol-
lowing.

u′1(x) = 1− x2

2
− x3

6
− x4

24
− x5

120
− x6

720
, (18)

u′2(x) = x+
x2

2
− x3

24
− x4

64
− x5

640
− x6

7680
− x7

107520
− x8

2580480
, (19)

u′3(x) = 2x+
3x2

2
+

x3

8
− 11x4

256
− 53x5

5120
− 23x6

40960
+

3x7

32768
+

575x8

22020096

+
297x9

73400320
+

28031x10

59454259200
+

71x11

1703116800
+

43x12

14863564800

+
x13

6038323200
+

x14

138726604800
+

x15

6242697216000
, (20)



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 7 of 16

Up to the n′th order problem, and by solving these problems and substituting them
in (??) , we have the HPM approximate solutions of order five

u5(x) = 1 + x+
x2

2
+

x3

6
+

x4

24
+

x5

120
− 42397x7

82575360
− 103239x8

1174405120
− 4027213x9

507343011840

− 16426909x10

81174881894400
+

1159751881x11

42860337640243200
+

3745787009x12

822918482692669440

+
29180882551x13

71319601833364684800
+

41682262283x14

1996948851334211174400

+
2848076657x15

17972539662007900569600
− 2620331837x16

41080090656018058444800

− 150638876327x17

26072164203019461092966400
− 9022206041257x18

28157937339261017980403712000

− 15360137046493x19

1070001618891918683255341056000
− 755701829501x20

1426668825189224911007121408000

− 101051x21

6290388881702765199360000
− 28139x23

15096933316086636478464000

− 19x24

3472294662699926390046720000
− 6079x22

166503640169427039682560000

− x25

1266875523028249214976000000
. (21)

This leads to u = ex. as limn→∞ ũn(x). Table 1 provides a comparative analysis of
the absolute errors obtained using the HPM on various orders: third-order, fifth-order,
and seventh-order. Absolute errors are assessed at different values of x, and presented in
graphical form in Fig 1. For example, at x = 0.2, the third-order HPM shows an error
of 7.40 × 10−8, while the fifth order significantly reduces it to 9.83 × 10−8, and the sev-
enth order further refines the result to an absolute error of 9.83 × 10−12. Similarly, for
larger values of x, the error decreases consistently as the order increases. The substantial



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 8 of 16

decrease in absolute error as the order increases highlights the effectiveness of HPM in
approximating solutions. These results indicate that, for practical applications that need
high precision, a seventh-order approximation or higher is recommended. The compar-
ison clearly shows that increasing the order of HPM leads to a significant reduction in
error. The third-order approximation provides a reasonable estimate, while the fifth-order
greatly improves accuracy, and the seventh-order results in an almost negligible error.
Therefore, for greater accuracy in solving non-linear problems, using a higher-order HPM
is advantageous. To improve the accuracy of the HPM approximation, we employ MHPM.
This method extends the original HPM to overcome its limitations, such as the number
of terms and computational complexity. The techniques we will consider include the Padé
approximation, Laplace transformation, and ultimately inverse Laplace, as follows.

L (ũ1(t)) = − 671

512s8
− 263

128s7
− 53

32s6
+

1

s4
+

1

s3
+

1

s2
+

1

s
, (22)

The use of s = 1
z , leads to

L (ũ1(t)) = z + z2 + z3 + z4 − 53z6

32
− 263z7

128
− 671z8

512
, (23)

The Pade approximation of order
[
3
3

]
in terms of x = 1

s , gives.[
3

3

]
=

1(
1− 1

s

)
s
, (24)

The modified approximation solution u1(x) = ex is obtained by applying the inverse
Laplace transform to the

[
3
3

]
Pade approximates.

7'th HPM absolute Error 

0.0 0.2 0.4 0.6 0.8 1.0

0

1×10
-6

2×10
-6

3×10
-6

4×10
-6

5×10
-6

5'th order HPM Absolute Error  

0.0 0.2 0.4 0.6 0.8 1.0

0.0000

0.0002

0.0004

0.0006

0.0008

0.0010

0.0012

Figure 1: Absolute error resulted from the HPM procedure for example 1

3.2. Example 2 The second example considered in this study is the fol-
lowing non-linear DVIDE [46]

.

u(x) = u
(x
2

)
− 3

2
sinx− x

2
− cos

(x
2

)
+

∫ x

0
u2

(s
2

)
ds, u(0) = 1. (25)



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 9 of 16

7'th HPM Absolute Error 

0.0 0.2 0.4 0.6 0.8 1.0

0

2×10
-10

4×10
-10

6×10
-10 5'th order HPM absolute error 

0.0 0.2 0.4 0.6 0.8 1.0

0

2×10
-9

4×10
-9

6×10
-9

8×10
-9

Figure 2: Absolute error resulted from the HPM procedure for example 2

Table 1: Numerical result of example 1

x 3th order HPM 5th order HPM 7th order HPM
Absolute Error Absolute Error Absolute Error

0.0 0.0000000000 0.0000000000 0.00
0.2 7.40× 10−8 9.83× 10−8 5.14× 10−12

0.4 1.31× 10−3 6.93× 10−6 2.75× 10−9

0.6 7.33× 10−3 8.67× 10−5 1.10× 10−7

0.8 2.55× 10−2 5.34× 10−4 1.53× 10−6

1.0 6.85× 10−2 2.22× 10−3 1.19× 10−5

Following the same process in example one, we have the 5th -order HPM approximate
solution.

L{u(x)} = u′(x),

N(u) = −u
(x
2

)
+

3

2

6∑
k=0

(−1)k
x2k+1

(2k + 1)!
+

x

2
+

6∑
k=0

(−1)k
(
x
2

)2k
(2k)!

−
∫ x

0
u2

(s
2

)
ds. (26)

Then we construct the homotopy equation.

u′(x)−u0(x)+p

[
−u

(x
2

)
+

3

2

6∑
k=0

(−1)k
x2k+1

(2k + 1)!
+

x

2
+

6∑
k=0

(−1)k
(
x
2

)2k
(2k)!

−
∫ x

0
u2

(s
2

)
ds

]
.

(27)
Using HPM, we expand u(x) as a perturbation series in p

u(x) = u0(x) + pu1(x) + p2u2(x) + p3u3(x) + p4u4(x) + p5u5(x) + · · · (28)

By substituting this series into the provided equation, we arrange the terms in accordance
with the powers of p, and we have: The zeroth order approximation that is obtained by
setting p = 0, is as follows

u′0(x) = 0, u0(0) = 1, (29)



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 10 of 16

which has the following solutions u0(x) = 1. Moreover, the first, second and other or-
ders of approximations will be obtained by substituting into the homotopy equations and
matching the terms of power p to get the following.

u′1(x) =
x

2
+
x2

8
− x4

384
+

x6

46080
− x8

10321920
+

x10

3715891200
−3

2

(
x− x3

6
+

x5

120
− x7

5040
+

x9

362880

)
,

(30)

u′2(x) = −x2

8
− 5x3

64
+

5x4

768
+

19x5

12288
− 7x6

184320
− 383x7

41287680
+

5x8

33030144

+
307x9

9512681472
− 97x10

237817036800
− 6143x11

83711596953600
+

x12

502269581721600
, (31)

u′3(x) = − x3

192
− 47x4

12288
+

329x5

122880
− 941x6

4718592
− 22217x7

165150720
+

471809x8

84557168640
+

1963291x9

761014517760

− 22998401x10

487049291366400
− 47100961x11

1785847401676800
+

37961641x12

164583696538533888
+

9442463747x13

53489701375023513600

− 25965661x14

34038900875014963200
− 1419107x15

1687625794584576000
+

24349x16

18001341808902144000

+
639917x17

229517108063502336000
− 17x18

11572291162865664000
− 2085431x19

337725745297071538176000

+
51622011x21

51843864409638174720000
+

x22

209034461299661120471040000

− 1

2281130034024079687680000
+

x23

161175523684005374412718080000
, (32)

Up to the n′th order problem, and by solving these problems and substituting them
in (??) , we have the HPM approximate solutions of order five

ũ(x) = 1− x2

2
+

x4

24
− x6

720
+

x7

5284823040
+

133153x8

5368709120
− 5419361x9

389639433093120

− 876413917x10

3117115464744960
+

6295085773x11

8777797148721807360
− 1814058456809x12

2106671315693233766400

− 130118761434241x13

7011002138627081974579200
+

455479147538587x14

130872039921038863525478400



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 11 of 16

Table 2: Numerical result of example 2

x 3th order HPM 5th order HPM 7th order HPM
Absolute Error Absolute Error Absolute Error

0.0 0.0000000000 0.0000000000 0.00
0.2 5.41× 10−9 5.11× 10−15 5.14× 10−12

0.4 1.52× 10−7 3.88× 10−12 5.26× 10−14

0.6 7.59× 10−7 1.71× 10−10 6.81× 10−12

0.8 5.213× 10−7 2.55× 10−9 2.15× 10−10

1.0 9.88× 10−6 2.15× 10−8 3.12× 10−9

+
19978465659899683x15

137058718171851609801228288000
− 716086078759136753x16

96489337592983533300064714752000
(33)

This converges to the exact solution ũ(x) = cosx, as in limn→∞ ũi(t). Table 2 provides
a comparative analysis of the absolute errors obtained using the HPM on various orders:
third-order, fifth-order, and seventh-order. Absolute errors are assessed at different val-
ues of x and plotted in Fig 2. For example, at x = 0.2, the third-order HPM shows an
error of 5.41 × 10−9, whereas the fifth order reduces it significantly to 5.11 × 10−15, and
the seventh order further refines the result to an absolute error of zero. Similarly, for
larger values of x, the error decreases consistently as the order increases. The substantial
decrease in absolute error as the order increases highlights the effectiveness of HPM in
approximating solutions. These results indicate that, for practical applications that need
high precision, a seventh-order approximation or higher is recommended. The compar-
ison clearly shows that increasing the order of HPM leads to a significant reduction in
error. The third-order approximation provides a reasonable estimate, while the fifth-order
greatly improves accuracy, and the seventh-order results in an almost negligible error.
Therefore, for greater accuracy in solving non-linear problems, using a higher-order HPM
is advantageous. Therefore, to enhance the precision of the HPM procedure, we will begin
by applying the Laplace transformation to the initial terms in the HPM series solutions.
Next, we will utilize the Pade approximants and, finally, we will conclude by implementing
the inverse Laplace transformation. The process is described below.

L (ũ1(t)) =
1

34359738368s10
+

1

s9
− 1

s7
+

1

s5
− 1

s3
+

1

s
, (34)

Use s = 1
z , leads to

L (ũ2(t)) = z − z3 + z5 − z7 + z9 +
z10

34359738368
, (35)



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 12 of 16

The Pade approximates of order
[
4
4

]
in term of x = 1

s , gives[
3

3

]
=

1(
1 + 1

s2

)
s
. (36)

The exact solutions ũ(x) = cosx are obtained by applying the inverse Laplace transform
to the

[
3
3

]
Pade approximate.

4. Conclusions

The present study holds considerable importance as it proposes a new and effective
analytical method for addressing Delay Volterra Integro-Differential Equations (DVIDEs),
which are critical for modeling systems influenced by memory and time-delay effects in
real-world scenarios. By enhancing the Homotopy Perturbation Method (HPM) and incor-
porating both the Laplace transformation and Padé approximants, the suggested approach
successfully overcomes common challenges associated with nonlinearity and limited con-
vergence in conventional methods. This combined technique not only streamlines the solu-
tion process but also improves accuracy, minimizes computational workload, and achieves
rapid convergence with fewer terms. Consequently, it provides a dependable and practical
solution for researchers and engineers dealing with mathematical models across various
domains such as biological systems, control theory, population models, and engineering
problems involving time delays and hereditary behavior. This work advances the field of
analytical techniques for integro-differential equations and opens new avenues for tackling
more complex or multidimensional problems.

References

[1] Vasile Marinca, Nicolae Herişanu, and Iacob Nemeş. Optimal homotopy asymptotic
method with application to thin film flow. Open Physics, 6(3):648–653, 2008.

[2] Vasile Marinca and Nicolae Herişanu. Application of optimal homotopy asymptotic
method for solving nonlinear equations arising in heat transfer. International com-
munications in heat and mass transfer, 35(6):710–715, 2008.

[3] Praveen Agarwal, Muhammad Akbar, Rashid Nawaz, and Mohamed Jleli. Solutions of
system of volterra integro-differential equations using optimal homotopy asymptotic
method. Mathematical Methods in the Applied Sciences, 44(3):2671–2681, 2021.

[4] N Ratib Anakira, AK Alomari, and Ishak Hashim. Application of optimal homotopy
asymptotic method for solving linear delay differential equations. In AIP Conference
Proceedings, volume 1571, pages 1013–1019. American Institute of Physics, 2013.

[5] N Ratib Anakira, AK Alomari, and Ishak Hashim. Numerical scheme for solv-
ing singular two-point boundary value problems. Journal of Applied Mathematics,
2013(1):468909, 2013.

[6] Remus-Daniel Ene and Nicolina Pop. Optimal homotopy asymptotic method for
an anharmonic oscillator: application to the chen system. Mathematics, 11(5):1124,
2023.



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 13 of 16

[7] Ali F Jameel, Akram H Shather, NR Anakira, AK Alomari, and Azizan Saaban.
Comparison for the approximate solution of the second-order fuzzy nonlinear differen-
tial equation with fuzzy initial conditions. Mathematics and Statistics, 8(5):527–534,
2020.

[8] Taiye Oyedepo Ayinde, Matthew Olanrewaju Oluwayemi, Muhammed Abdullahi,
Johnson Adekunle Osilagun, and Lukman Olalekan Ahmed. Homotopy perturba-
tion technique for fractional volterra and fredholm integro differential equations. In
2023 International Conference on Science, Engineering and Business for Sustainable
Development Goals (SEB-SDG), volume 1, pages 1–6. IEEE, 2023.

[9] Samad Noeiaghdam, Aliona Dreglea, Jihuan He, Zakieh Avazzadeh, Muhammad Sule-
man, Mohammad Ali Fariborzi Araghi, Denis N Sidorov, and Nikolai Sidorov. Error
estimation of the homotopy perturbation method to solve second kind volterra in-
tegral equations with piecewise smooth kernels: Application of the cadna library.
Symmetry, 12(10):1730, 2020.

[10] Moustafa El-Shahed. Application of he’s homotopy perturbation method to volterra’s
integro-differential equation. International Journal of Nonlinear Sciences and Numer-
ical Simulation, 6(2):163–168, 2005.

[11] Nidal Anakira, Adel Almalki, MJ Mohammed, Safwat Hamad, Osoma Oqilat, and
Ala Amourah. Analytical approaches for computing exact solutions to system of
volterra integro-differential equations. WSEAS Trans. Math, 23:400–407, 2024.

[12] S Alao, FS Akinboro, FO Akinpelu, and R Oderinu. Numerical solution of integro-
differential equation using adomian decomposition and variational iteration methods.
IOSR Journal of Mathematics, 10(4):18–22, 2014.

[13] Mehdi Dehghan, Mohammad Shakourifar, and Asgar Hamidi. The solution of linear
and nonlinear systems of volterra functional equations using adomian–pade technique.
Chaos, Solitons & Fractals, 39(5):2509–2521, 2009.

[14] Nidal Anakira, Gada Bani-Hani, Osama Ababneh, Ali Jameel, and Khamis Al-
Kalbani. Modified adomian decomposition method for solving volterra integro-
differential equations. In The International Arab Conference on Mathematics and
Computations, pages 335–341. Springer, 2023.

[15] MA Fariborzi Araghi and Sh Sadigh Behzadi. Solving nonlinear volterra— fredholm
integro-differential equations using the modified adomian decomposition method.
Computational Methods in Applied Mathematics, 9(4):321–331, 2009.

[16] Xindong Zhang, Bo Tang, and Yinnian He. Homotopy analysis method for higher-
order fractional integro-differential equations. Computers & mathematics with appli-
cations, 62(8):3194–3203, 2011.

[17] AF Jameel, Azizan Saaban, SA Altaie, NR Anakira, AK Alomari, and N Ah-
mad. Solving first order nonlinear fuzzy differential equations using optimal homo-
topy asymptotic method. International Journal of Pure and Applied Mathematics,
118(1):49–64, 2018.

[18] Ali F Jameel, Akram H Shather, NR Anakira, AK Alomari, and Azizan Saaban.
Comparison for the approximate solution of the second-order fuzzy nonlinear differen-
tial equation with fuzzy initial conditions. Mathematics and Statistics, 8(5):527–534,



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 14 of 16

2020.
[19] Nidal Ratib Anakira, AF Jameel, Abedel-Karrem Alomari, Azizan Saaban, Moham-

mad Almahameed, and IHM Hashim. Approximate solutions of multi-pantograph
type delay differential equations using multistage optimal homotopy asymptotic
method. Journal of Mathematical and Fundamental Sciences, 2018.

[20] FM Alharbi. Numerical solutions of an integro-differential equation with smooth and
singular kernels. International Journal of Mathematical Analysis, 13(12):573–586,
2019.

[21] Ivo Babuška, Fabio Nobile, and Raúl Tempone. A stochastic collocation method
for elliptic partial differential equations with random input data. SIAM Journal on
Numerical Analysis, 45(3):1005–1034, 2007.

[22] H Ghaneai, MM Hosseini, and Syed Tauseef Mohyud-Din. Modified variational it-
eration method for solving a neutral functional-differential equation with propor-
tional delays. International Journal of Numerical Methods for Heat & Fluid Flow,
22(8):1086–1095, 2012.

[23] TL Yookesh, ED Boobalan, and TP Latchoumi. Variational iteration method to
deal with time delay differential equations under uncertainty conditions. In 2020
International Conference on Emerging Smart Computing and Informatics (ESCI),
pages 252–256. IEEE, 2020.

[24] Ahmet Yıldırım. Applying heś variational iteration method for solving differential-
difference equation. Mathematical Problems in Engineering, 2008(1):869614, 2008.

[25] Ali Fareed Jameel, Nidal Ratib Anakira, AK Alomari, M Al-Mahameed, and Azizan
Saaban. A new approximate solution of the fuzzy delay differential equations. Interna-
tional Journal of Mathematical Modelling and Numerical Optimisation, 9(3):221–240,
2019.

[26] Brijendra Kumar Chaurasiya and Avadhesh Kumar. Analysis of positive solutions for
the fractional derivative with delay and integral boundary conditions. Gulf Journal
of Mathematics, 19(2):315–327, 2025.

[27] Zhor Mellah, El Bekkaye Mermri, and Mohammed Bouchlaghem. A numerical ap-
proach for solving a semilinear obstacle problem on the boundary. Gulf Journal of
Mathematics, 19(2):20–35, 2025.

[28] Mouhssine Zakaria and Abdelaziz Moujahid. Computational spectral method for solv-
ing two-dimensional riesz multi-term time-fractional diffusion equation. Gulf Journal
of Mathematics, 19(2):181–199, 2025.

[29] Ali Jameel, NR Anakira, AK Alomari, Ishak Hashim, and MA Shakhatreh. Numerical
solution of n’th order fuzzy initial value problems by six stages. Journal of nonlinear
science and applications, 2016.

[30] N Rabit Anakira, AK Alomari, AF Jameel, and Ishak Hashim. Multistage optimal
homotopy asymptotic method for solving initial-value problems. J. Nonlinear Sci.
Appl, 9(4):1826–1843, 2016.

[31] AF Jameel, NR Anakira, MM Rashidi, AK Alomari, A Saaban, and MA Shakhatreh.
Differential transformation method for solving high order fuzzy initial value problems.
Italian Journal of Pure and Applied Mathematics, 39:194–208, 2018.



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 15 of 16

[32] Erkan Cimen and Sabahattin Yatar. Numerical solution of volterra integro-differential
equation with delay. J. Math. Comput. Sci, 20(3):255–263, 2020.

[33] Rohul Amin, Kamal Shah, Muhammad Asif, and Imran Khan. Efficient numerical
technique for solution of delay volterra-fredholm integral equations using haar wavelet.
Heliyon, 6(10), 2020.

[34] Erkan Cimen and Sabahattin Yatar. Numerical solution of volterra integro-differential
equation with delay. J. Math. Comput. Sci, 20(3):255–263, 2020.

[35] Emiko Ishiwata and Yoshiaki Muroya. On collocation methods for delay differential
and volterra integral equations with proportional delay. Frontiers of mathematics in
China, 4(1):89–111, 2009.

[36] Aidouni Yamina. Numerical Solution of Integro-Delay Differential Equations on a
Half Line. PhD thesis, University BBA, 2023.

[37] Ghada H Ibraheem, Mustafa Turkyilmazoglu, and MA Al-Jawary. Novel approxi-
mate solution for fractional differential equations by the optimal variational iteration
method. Journal of Computational Science, 64:101841, 2022.

[38] Adedapo Ismaila Alaje, Morufu Oyedunsi Olayiwola, Kamilu Adewale Adedokun,
Joseph Adeleke Adedeji, and Asimiyu Olamilekan Oladapo. Modified homotopy
perturbation method and its application to analytical solitons of fractional-order
korteweg–de vries equation. Beni-Suef University Journal of Basic and Applied Sci-
ences, 11(1):139, 2022.

[39] Adedapo Ismaila Alaje and Morufu Oyedunsi Olayiwola. A fractional-order mathe-
matical model for examining the spatiotemporal spread of covid-19 in the presence
of vaccine distribution. Healthcare Analytics, 4:100230, 2023.

[40] Morufu Oyedunsi Olayiwola and Adedapo Ismaila Alaje. Mathematical analysis of
intrahost spread and control of dengue virus: unraveling the crucial role of antigenic
immunity. Franklin Open, 7:100117, 2024.

[41] Mutairu Kayode Kolawole, Morufu Oyedunsi Olayiwola, Adedapo Ismaila Alaje,
Hammed Ololade Adekunle, and Kazeem Abidoye Odeyemi. Conceptual analysis
of the combined effects of vaccination, therapeutic actions, and human subjection to
physical constraint in reducing the prevalence of covid-19 using the homotopy pertur-
bation method. Beni-Suef University Journal of Basic and Applied Sciences, 12(1):10,
2023.

[42] M34108671396 Turkyilmazoglu. Is homotopy perturbation method the traditional
taylor series expansion. Hacettepe journal of mathematics and statistics, 44(3):651–
657, 2015.

[43] Mustafa Turkyilmazoglu. Optimization by the convergence control parameter in iter-
ative methods. Journal of Applied Mathematics and Computational Mechanics, 23(2),
2024.

[44] Shadi Al-Ahmad, Mustafa Mamat, Nidal Anakira, and Rami Alahmad. Modified dif-
ferential transformation method for solving classes of non-linear differential equations.
TWMS Journal of Applied and Engineering Mathematics, 2022.

[45] Shadi Al-Ahmad, Nidal Ratib Anakira, Mustafa Mamat, Ali Fareed Jameel, Rami
Alahmad, and Abedel-Karrem Alomari. Accurate approximate solution of classes



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6462 16 of 16

of boundary value problems using modified differential transform method. TWMS
Journal of Applied and Engineering Mathematics, 2022.

[46] Fazli Hadi, Rohul Amin, Ilyas Khan, J Alzahrani, KS Nisar, Amnah S Al-Johani,
and Elsayed Tag Eldin. Numerical solutions of nonlinear delay integro-differential
equations using haar wavelet collocation method. Fractals, 31(02):2340039, 2023.