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

A Numerical Comparison Between the Standard and
Modified Versions of OHAM Procedure for Solving

Volterra Integro-Delay Differential Equations
Nidal Anakira1,∗, Sameer Bawaneh2, Areen Al-khateeb3, Ala Amourah1,
Tala Sasa4

1 Mathematics Education Program, Faculty of Education and Arts, Sohar University,
Sohar 311, Oman
2 Department of Computer Science, Faculty of Information Technology, Jadara University,
Irbid, 21110, Jordan
3 Department of Mathematics, Faculty of Science and Technology, Jadara University,
Irbid, 21110, Jordan
4 Department of Mathematics, Faculty of Science, Private Applied Science University,
Amman, Jordan

Abstract. In this paper, we formulate and enhance a robust semi-analytical method, known as the
Optimal Homotopy Asymptotic Method (OHAM), for solving Volterra delay integro-differential
equations (VDIDEs). A comparative analysis between the standard OHAM and its modified ver-
sion is presented, emphasizing the improvements introduced through the modification. The modi-
fication is based on a refined construction of the auxiliary function H(p), which plays a crucial role
in enhancing both the accuracy and the convergence of the method. The effectiveness and efficiency
of the proposed technique are demonstrated through the solution of various numerical problems.
Notably, the modified OHAM achieves higher accuracy within a single order of iteration, in con-
trast to the four orders required by the standard OHAM. This advancement significantly reduces
computational effort, simplifies calculations, and decreases overall time consumption, thereby es-
tablishing the modified OHAM as a more efficient and practical tool for addressing such equations.
2020 Mathematics Subject Classifications: 26A33, 34A08
Key Words and Phrases: OHAM, MOHAM, nonlinear delay Volterra integro-differential equa-
tions

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.7040

Email addresses: nanakira@su.edu.om (N. Anakira),
s.bawaneh@jadara.edu.jo (S. Bawaneh), a.k@jadara.edu.jo (A. Al-khateeb),
aamourah@su.edu.om (A. Amourah), 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 (4) (2025), 7040 2 of 16

1. Introduction

Differential and integral equations play a central role in science and engineering, with
applications spanning control theory, economics, electrical engineering, medicine, and
many other disciplines. Since exact solutions to nonlinear systems are often difficult
to obtain, numerous researchers have contributed alternative approaches for construct-
ing approximate solutions. For example, Yang et al. [1] applied the modified reproducing
kernel method to linear Volterra integral equations, while Anakira et al. [2] and Raftari [3]
successfully implemented the homotopy perturbation method (HPM) for Volterra integro-
differential and delay equations. The Chebyshev wavelet method was developed and ap-
plied by Biazar and Ebrahimi [4] and later extended by Sahu and Ray [5] for Lane–Emden
type equations, Anakira et al. [6] proposed a modified version for handling nonlinear
delay integro-differential problems. whereas Dawar et al. [7] further advanced the resid-
ual power series method with an improved formulation In the context of homotopy-based
strategies, Hashim et al. [8], Jameel et al. [9, 10], and Nawaz et al. [11] demonstrated the
flexibility of the OHAM for fractional, fuzzy, and integro-differential equations. Jameel
et al. [12] applied the differential transformation method (DTM) for high-order fuzzy
problems, while Matinfar et al. [13] employed the homotopy analysis method (HAM) to
systems of Volterra integral equations. Hamoud et al. [14] investigated the modified vari-
ational iteration method, and Hong et al. [15] applied the relaxed Monte Carlo method
for Volterra integral systems. Other notable contributions include Agarwal et al. [16],
who explored collocation methods for stochastic equations, while Anakira et al. [17, 18]
introduced Multistage MOHAM for solving fuzzy and delay differential equations. Finally,
Jameel et al. [19] developed a six-stage Runge–Kutta method of order five for n′th order
fuzzy initial value problems. Together, these studies underscore the rich diversity of tech-
niques proposed by different authors, each of whom has advanced the state of the art in
approximating solutions for complex differential and integro-differential equations [20–24].

This research employs the OHAM [25, 26] to solve the Volterra delay integro-differential
equations (VDIDE). The motivation stems from the need to achieve highly accurate ap-
proximate solutions for such systems by employing a robust semi-analytical methodology.
A key advantage of OHAM lies in its flexible convergence, which provides an effective
framework for constructing and controlling approximation series. employed it to solve
Klein–Gordon equations [27], Sheikholeslami et al. used it to investigate laminar viscous
flow [28]; Hashmi et al. applied it to Fredholm integral equations [29], and Nawaz et
al. derived optimal solutions for fractional Zakharov–Kuznetsov equations [30], three-
dimensional Volterra integral equations [31], and fractional integro-differential equations
[32].

In this study, we build on these foundations by modifying the standard OHAM to
overcome some of its inherent difficulties and limitations. The proposed modification
enhances convergence and accuracy while reducing computational effort. We then perform
a comparative study between the standard OHAM and its modified version, supported
by two numerical examples, to clearly demonstrate the efficiency and reliability of the
improved approach in solving systems of VIEs.



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

Modified OHAM for Volterra Integro-Delay Differential Equations

We consider the general nonlinear Volterra integro-delay differential equation of the
form

L[u(x)] + g(x) +N [u(x)] +

∫ x

0
K(x, s)u(s− τ) ds = 0, x ∈ [0, T ], τ > 0, (1)

subject to the boundary/initial conditions

u(0) = u0, B

(
du

dx

)
x=0

= 0, (2)

where L is a linear operator, N is a nonlinear operator, K(x, s) is a known kernel, τ is a
constant time-delay, and g(x) is a given source function.

According to the basic idea of OHAM, we construct the homotopy

(1− p)L[v(x, p)] = H(p)
[
L(v(x, p)) + g(x) +N(v(x, p)) +

∫ x

0
K(x, s) v(s− τ, p) ds

]
, (3)

with the boundary operator;

u(0) = u0, B

(
du

dx

)
x=0

= 0, (4)

where p ∈ [0, 1] is the embedding parameter, v(x, p) is the unknown function, and H(p) is
a non-zero auxiliary function with H(0) = 0, H(1) = 1.

To improve convergence, we introduce a modified auxiliary function in the form

H(q) =

N∑
n=1

(
M∑
k=0

cnkx
k

)
qn, (5)

where cnk are the convergence control parameters to be determined.

Expanding v(x, p) into a power series of p, we have

v(x, p) = u0(x) +

∞∑
k=1

uk(x,C1, C2, . . . , Cm) pk. (6)

When p→ 1, the approximate solution becomes

u(x) ≈ u0(x) +
m∑
k=1

uk(x,C1, C2, . . . , Cm). (7)

The residual for the truncated m-term approximation is given by

R(x;C1, C2, . . . , Cm) = L[ũ(x)] + g(x) +N [ũ(x)] +

∫ x

0
K(x, s) ũ(s− τ) ds, (8)



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

where

ũ(x) = u0(x) +

m∑
k=1

uk(x,C1, C2, . . . , Cm). (9)

The optimal values of the unknown convergence control constants C1, C2, . . . , Cm are ob-
tained by minimizing the functional value.

J(C1, C2, . . . , Cm) =

∫ b

a
R2(x;C1, C2, . . . , Cm) dx, (10)

which leads to the conditions

∂J

∂C1
=

∂J

∂C2
= · · · = ∂J

∂Cm
= 0. (11)

Thus, with the modified auxiliary function, the OHAM procedure provides an approx-
imate analytical solution to Volterra integro-delay differential equations with improved
convergence properties.

2. Applications

In this section, we present two examples with known exact solutions to evaluate the
performance and accuracy of the MOHAM algorithm in comparison with the standard
OHAM.
Example 1 As a first step, we consider the following nonlinear VDIDE as a test case.
It is solved using OHAM and then with MOHAM, enabling a direct comparison that
demonstrates the efficiency, accuracy, and reliability of the modified approach. [5]

u′(x) = u2
(x
2

)
− ex +

∫ x

0
u2(

t

2
)dt+ 1, u(x) = ex, x ≤ 0, (12)

with the exact solution

u(x) = ex. (13)

2.1. OHAM solution

Following the OHAM formulation described in [25], we begin by introducing the cor-
responding operators that will be used in constructing the solution.

L[v(x, p)] =
dv(x, p)

dx
= 1,

N [v(x, p)] = u′(x)− u2
(x
2

)
− ψ(x) +

∫ x

0
u2(

t

2
)dt+ 1, v(0, p) = 1,

(14)



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

where ψ(x) is the expansion Taylor series of ex with respect to p, which can be written as

ψ(x; p) =
∞∑
k=0

( x
k!

)k
pk.

Then, we construct a family of homtopy equations h(u(x; p) : R− > [0, 1] as follows

(1− p)L[v(x, p)] = H(p)
[
L(v(x, p)) + g(x) +N(v(x, p))

]
, (15)

substituting L[u(x)], we have

u′(x; p)− u0(x) +
m∑
i=1

Cmq
m

(
(u′(x)− u2

(x
2

)
− ψ(x) +

∫ x

0
u2(

t

2
)dt+ 1

)
= 0. (16)

Thus, as p varies from 0 to 1, the solution approach from u0(x) to u(x), where u0(x) is
the zeroth -order problem that can be obtained from the solution of initial guess

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

Then, by expanding u(x; p; ci), i = 1, ..., 4 in Taylor series about p, we have

u(x; p; ci) = u0(x) +

m∑
k=1

uk(x,C1, C2, . . . , Cm) pk. (18)

Substituting Eq.(18) into Eq. (16) and by equating the coefficient of like power p, yields
The first-order deformation problem

u′1(x,C1) = c1(−x)−
1

4

(
c1x

2
)
+
c1x

3

12
+
c1x

4

24
+
c1x

5

120
+
c1x

6

720
+
c1x

7

5040
, u0(0) = 0, (19)

The second-order deformation problem is

u′2(x,C1, C2) = c1(−x)−
1

4

(
c1x

2
)
+
c1x

3

12
+
c1x

4

24
+
c1x

5

120
+
c1x

6

720
+
c1x

6

720

+
c1x

7

5040
c21(−x) +

5

16
c21x

3 +
11

128
c21x

4 +
31c21x

5

3840
+
c21x

6

1280
+
c21x

7

7680

− 5c21x
8

1032192
− c21x

9

3440640
− c21x

10

103219200
− c2x− 1

4

(
c2x

2
)
+
c2x

3

12

+
c2x

4

24
+
c2x

5

120
+
c2x

6

720
+
c2x

7

5040
, u2(0) = 0. (20)

The third-order deformation problem

u′3(x,C1, C2, C3) = − 5c21x
8

172032
+
c21x

7

1280
+

3

640
c21x

6 +
31

640
c21x

5 +
33

64
c21x

4 +
15

8
c21x

3 − 6c21x

− c21x
10

17203200
− c21x

9

573440
+

141

256
c31x

4 +
25

8
c31x

3 +
3

2
c31x

2 − 6c31x



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

8509c31x
9

1981808640
+

283c31x
8

3932160
+

509c31x
7

573440
− 743c31x

6

122880
− 171c31x

5

5120

− 13801c31x
12

871995801600
− 12659c31x

11

108999475200
− 3509c31x

10

8808038400
+ . . . , u3(0) = 0,

(21)

and the fourth order deformation problem

u′4(x,C1, C2, C3, C4) =
c1x

7

5040
+
c1x

6

720
+
c1x

5

120
+
c1x

4

24
+
c1x

3

12
− c1x

2

4
− 3c21x− c1x

− 5c21x
8

344064
+
c21x

7

2560
+

3c21x
6

1280
+

31c21x
5

1280
+

33

128
c21x

4 +
15

16
c21x

3

− c21x
10

34406400
− c21x

9

1146880
+

25

16
c31x

3 +
3

4
c31x

2 − 3c31x

283c31x
8

7864320
+

509c31x
7

1146880
− 743c31x

6

245760
− 171c31x

5

10240
+

141

512
c31x

4

− 12659c31x
11

217998950400
− 3509c31x

10

17616076800
+

8509c31x
9

3963617280
+ . . . , u4(0) = 0

(22)

By solving Eqs. (17), (19), (20), (21) and (22) and substituting them into Eq.(18), the
fourt-order approximate solution by OHAM for p = 1, is:

ũ(x,C1, C2, C3, C4) = u0(x) + u1(x,C1) + u2(x,C1, C2)

+u3(x,C1, C2, C3) + +u4(x,C1, C2, C3, c4). (23)

By using the proposed method of section 2 on [0, 1], we use the residual error,

R = ũ′ (x,C1, · · · , C4)− ũ2
(x
2
, C1, · · · , C4

)
+ex −

∫ x

0
u2(

t

2
, C1, · · · , C4)dt− 1. (24)

The Less Square error can be formed as

J(C1, · · · , C4) =

∫ 1

0
R2dx, (25)

and

∂J(C1, · · · , C4)

∂C1
=
∂J(C1, C2, C3)

∂C2
=
∂J(C1, C2, C3)

∂C3
=
∂J(C1, C2, C3, C4)

∂C4
= 0. (26)

Thus, the following optimal values of Ci’s are obtained:

C1 = −1.03257, C2 = 3.75206× 10−6, C3 = −8.89991× 10−7 C4 = −4.48187× 10−6.



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

In this case, our approximate solution is

ũ(x,C1, · · · , C4) = 1 + x+ 0.5000016432167662x2 + 0.16667895652842996x3

58556x4 + 0.00830125122279865x5 + 0.0014533846838279182x6

+0.000205626800537755x7 + 0.000006767376975033287x8

−0.000001498570429380608x9 − 6.299736808070288× 10−7x10

−6.739802323017767× 10−8x11− 3.758169088190813× 10−9x12

−7.089767587924317× 10−11x13 + 9.59982558191048× 10−12x14

+1.309622277385672× 10−12x15 + 8.75719622963593× 10−14x16

+3.636632322436979× 10−15x17 + 9.172211052795001× 10−17x18

+1.443596286575417× 10−18x19 + 2.522219242033677× 10−20x20

+2.258535766371391× 10−22x21 (27)

As we observed, the accuracy of the OHAM strongly depends on the number of approxi-
mation terms used in the series solution. While increasing the number of terms generally
improves the precision of the obtained results and reduces the error, it also leads to heav-
ier computational work and longer calculation time, which can limit the efficiency of the
method when applied to complex nonlinear problems. These practical difficulties—namely
the high computational cost, extended time requirements, and the necessity of employing
higher-order approximation terms to achieve remarkably low errors motivated the modifi-
cation of OHAM. The proposed modification is built upon the construction of the auxiliary
function H(p) (Eq.(5)) in the standard OHAM. This formulation offers enhanced flexibility
in controlling the convergence of the series solution and enables the method to effectively
applied to a wider class of nonlinear problems with only one order of approximation. This
feature represents the main advantage of the proposed modification, as it significantly save
time, reduces effort, and minimizes computational calculations as we will see on next part.

2.2. MOHAM Solution

The MOHAM solution is obtained by constructing a family of homotopy equations
similar to the standard OHAM framework, but with an important modification in the de-
sign of the auxiliary function. This new construction of the auxiliary function introduces
additional flexibility through adjustable convergence-control parameters, which allows the
method to achieve accurate results using only a single term of the OHAM order of ap-
proximation. Based on Eq. (16), and by employing Eq. (5) as follows

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

N∑
n=1

(
M∑
k=0

cnkx
k

)
qn)

(
(u′(x)− u2

(x
2

)
− ψ(x) +

∫ x

0
u2(

t

2
)dt+ 1

)
= 0.

(28)

Thus, as p varies from 0 to 1, the solution approach from u0(x) to u(x), where u0(x) is
the zeroth -order problem that can be obtained from the solution of initial guess

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



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

Then, by expanding u(x; p; ci), i = 1, ..., 5, and m = 1 in Taylor series about p, we have

u(x; p; ci) = u0(x) +
m∑
k=1

uk(x,C1, C2, . . . , Ci) p
k. (30)

Substituting Eq.(30) into Eq. (28) and by equating the coefficient of like power p, yields
to the first order deformation problem

u′1(x) =
c5x

12

5040
+
c4x

11

5040
+
c5x

11

720
+
c3x

10

5040
+
c4x

10

720
+
c5x

10

120
+
c2x

9

5040
+
c3x

9

720
+
c4x

9

120

+
c5x

9

24
+
c1x

8

5040
+
c2x

8

720
+
c3x

8

120
+
c4x

8

24
+
c5x

8

12
+
c0x

7

5040
+
c1x

7

720
+
c2x

7

120
+
c3x

7

24

+
c4x

7

12
− c5x

7

4
+
c0x

6

720
+
c1x

6

120
+
c2x

6

24
+
c3x

6

12
− c5x

6 − c4x
6

4
+
c0x

5

120
+
c1x

5

24

+
c2x

5

12
− c4x

5 − c3x
5

4
+
c0x

4

24
+
c1x

4

12
− c3x

4 − c2x
4

4
+
c0x

3

12
− c2x

3 − c1x
3

4

−c1x2 −
c0x

2

4
− c0x (31)

By solving Eqs. (31) and (30), and substituting them into Eq.(30), we obtain MOHAM
approximate solution of order one as follows:

u(x, c1, . . . , c5) =
c5x

13

65520
+
c4x

12

60480
+
c5x

12

8640
+
c3x

11

55440
+
c4x

11

7920
+
c5x

11

1320
+
c2x

10

50400
+
c3x

10

7200

+
c4x

10

1200
+
c5x

10

240
+

c1x
9

45360
+
c2x

9

6480
+
c3x

9

1080
+
c4x

9

216
+
c5x

9

108
+

c0x
8

40320

+
c1x

8

5760
+
c2x

8

960
+
c3x

8

192
+
c4x

8

96
− c5x

8

32
+
c0x

7

5040
+
c1x

7

840
+
c2x

7

168
+
c3x

7

84

−c5x
7

7
− c4x

7

28
+
c0x

6

720
+
c1x

6

144
+
c2x

6

72
− c4x

6

6
− c3x

6

24
+
c0x

5

120
+
c1x

5

60

−c3x
5

5
− c2x

5

20
+
c0x

4

48
− c2x

4

4
− c1x

4

16
− c1x

3

3
− c0x

3

12
− c0x

2

2
+ x+ 1.

(32)

the residual error,

R = ũ′ (x,C1, · · · , C5)− ũ2
(x
2
, C1, · · · , C5

)
+ex −

∫ x

0
u2(

t

2
, C1, · · · , C5)dt− 1. (33)

The Less Square error can be formed as

J(C1, · · · , C5) =

∫ 1

0
R2dx, (34)



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

Table 1: Numerical result of example 1

x 4’th order OHAM First-order MOHAM
Absolute Error Absolute Error

0.0 0.0000000000 0.0000000000
0.2 5.41× 10−9 2.51× 10−7

0.4 1.52× 10−7 4.58× 10−7

0.6 7.59× 10−7 5.86× 10−7

0.8 5.213× 10−7 4.30× 10−7

1.0 9.88× 10−6 1.71× 10−6

and

∂J(C1, · · · , C5)

∂C1
= · · · = ∂J(C1, C2, C3, C5)

∂C5
= 0. (35)

Thus, the following optimal values of Ci’s are obtained:

C1 = −0.999749, C2 = −0.254198, C3 = −0.122836 C4 = 0.0566211 C5 = −0.0578005.

In this case, our approximate solution is

u(x) = −8.821804171278825× 10−7x13 − 5.753673082235477×1 0−6x12

−0.0000388548x11 − 0.000213946x10 − 0.00041755x9 + 0.00151758x8

+0.00330149x7 − 0.00973638x6 + 0.0201492x5 + 0.035809x4

+0.168045x3 + 0.499875x2 + x+ 1 (36)

The numerical results of Example 1 highlight the superior efficiency of MOHAM compared
with OHAM. As shown in Table 1, the fourth-order OHAM solution produces extremely
small errors (as low as 5.41 × 10−9 at x = 0.2), while the first-order MOHAM approxi-
mation already achieves comparable accuracy, with absolute errors ranging from 10−7 to
10−6. This indicates that MOHAM can provide reliable results with only a single term,
thereby reducing computational complexity and effort. Furthermore, the error distribution
presented in Figure 1 demonstrates that MOHAM ensures smooth and stable convergence
across the interval, whereas OHAM requires higher-order expansions to attain a similar
level of accuracy. Overall, these findings clearly confirm the advantage of MOHAM in
achieving high accuracy at lower approximation orders.

2.3. Example 2

The second example considered in this study is the following non-linear DVIDE

u(x) = u
(x
2

)
− 3

2
sinx− x

2
− cos

(x
2

)
+

∫ x

0
u2
(s
2

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



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

a) 0.0 0.2 0.4 0.6 0.8 1.0

0

5.0× 10
-7

1.0× 10
-6

1.5× 10
-6

b) 0.0 0.2 0.4 0.6 0.8 1.0

0

5.0× 10
-7

1.0× 10
-6

1.5× 10
-6

Figure 1: (a) OHAM absolute error resulted from fourth -order of approximation. (b) MOHAM absolute error
resulted from asingle -order of approximation foe example 1

with exact solution

u(x) = cos(x). (38)

According to the method which was described in above section, and by following the
same procedure that was applied on the first example, we have the fourth order OHAM
approximate solutions

ũ(x,C1, C2, C3) = 6.874157120897826× 10−37x26 − 8.463767623593088× 10−32x25

+2.184023696562375× 10−27x24 + 1.0133589762100683× 10−26x23

+3.460786587767694× 10−22x22 − 2.7806197391988356× 10−23x21

−2.9629995732517603× 10−19x20 + 5.564560645958861× 10−20x19

+1.3918670220163003× 10−16x18 − 9.700555588641355× 10−17x17



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

Table 2: Numerical result of example 2

x 4’th order OHAM First-order MOHAM
Absolute Error Absolute Error

0.0 0.00 0.00
0.2 6.06× 10−8 3.58× 10−7

0.4 1.81× 10−7 7.29× 10−7

0.6 2.36× 10−7 6.74× 10−7

0.8 1.47× 10−7 1.01× 10−6

1.0 2.57× 10−7 1.25× 10−6

−4.222517784937938× 10−14x16 + 4.202560325267928× 10−14x15

+9.028868136417469× 10−12x14 − 1.2725695013294234× 10−11x13

−1.3653035357656884× 10−9x12 + 2.555212476061491× 10−9x11

−2.6421049209252045× 10−7x10 − 2.710683819309949× 10−7x9

+0.0000239393x8 + 6.931541036055155× 10−8x7 − 0.00138905x6

+3.5759065280388524× 10−7x5 + 0.041668x4

+1.0621242323427538× 10−6x3 − 0.500002x2 + 1

ũ(x,C1, C2, C3) = 1− 0.499725609x2 − 0.003203684x3

+0.056106703x4 − 0.03141667x5

+0.031943943x6 − 0.0128500196x7

−0.0035950170x8 + 0.002988065x9

+0.000183446x10 − 0.000126185x11

−4.018124199× 10−6x12 + 2.579615120× 10
−6x13

+4.9331251480× 10−8x14 − 3.177894121× 10
−8x15

+1.93944674× 10−12x16.



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

a) 0.0 0.2 0.4 0.6 0.8 1.0

0

5.0× 10
-8

1.0× 10
-7

1.5× 10
-7

2.0× 10
-7

2.5× 10
-7

b) 0.0 0.2 0.4 0.6 0.8 1.0

0

5.0× 10
-7

1.0× 10
-6

1.5× 10
-6

Figure 2: (a) OHAM absolute error resulted from fourth -order of approximation. (b) MOHAM absolute error
resulted from asingle -order of approximation foe example 2

The results of Example 2 further confirm the robustness of MOHAM in comparison
with OHAM. As shown in Table 2, the fourth–order OHAM approximation attains errors
in the order of 10−8 to 10−7, while the first–order MOHAM already achieves comparable
accuracy, with errors ranging from 10−7 to 10−6. Although OHAM demonstrates slightly
higher precision at certain points (e.g., 6.06×10−8 at x = 0.2), MOHAM reaches this level
of accuracy with far fewer terms, thus reducing the computational effort. The error profiles
illustrated in Figure 2 indicate that OHAM’s residual exhibits oscillatory behavior, while
MOHAM maintains a smoother distribution of error across the interval. This comparison
highlights the advantage of MOHAM in producing reliable results at lower approximation
orders, striking a balance between computational efficiency and accuracy.

The Modified Optimal Homotopy Asymptotic Method (MOHAM) offers remarkable
advantages over the standard OHAM, especially in terms of the order of approximation.



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

Unlike the standard OHAM, which generally requires more than two approximation terms
to achieve acceptable accuracy, MOHAM attains highly accurate solutions with only a
single approximation term. This efficiency not only reduces the analytical effort but also
saves considerable time in the solution process. In addition, MOHAM minimizes com-
putational calculations, making it more practical for solving complex nonlinear problems.
These benefits are clearly demonstrated in the numerical results displayed in Tables 1 and
2, where MOHAM consistently provides accurate solutions with fewer terms and reduced
computational effort compared to the standard OHAM.

3. Conclusions

In this research, we conducted a comparative analysis between the conventional OHAM
and a newly suggested modification, which was applied for the first time to derive an
analytic approximate solution for VDIDEs. The modified approach introduced a new
series-based formulation for the auxiliary function H(q), specifically designed to improve
convergence, simplify implementation, and reduce computational complexity and cost,
while ensuring high accuracy and reliability. A key advantage of the modified technique
is its ability to deliver accurate results using only one order of approximation, in contrast
to the standard OHAM, which typically requires several terms to achieve similar levels of
accuracy. This highlights the improved efficiency, robustness, and practical applicability
of the modified version in addressing this class of VDIDEs.

Acknowledgements

The authors gratefully acknowledge Jadara University for supporting this research and
for encouraging scientific work.

References

[1] Li-Hong Yang, Hong-Ying Li, and Jing-Ran Wang. Solving a system of linear volterra
integral equations using the modified reproducing kernel method. 2013(1):196308,
2013.

[2] Behrouz Raftari. Numerical solutions of the linear volterra integro-differential equa-
tions: Homotopy perturbation method and finite difference method. World Applied
Sciences Journal, 9(07-12), 2010.

[3] Nidal Anakira, Ala Amourah, Adel Almalki, Abdullah Alsoboh, and Tala Sasa. Ex-
act solutions of nonlinear delay voltera integro-differential equations using modified
homotopy perturbation method. European Journal of Pure and Applied Mathematics,
18(3):6462–6462, 2025.

[4] Jafar Biazar and H Ebrahimi. Chebyshev wavelets approach for nonlinear systems of
volterra integral equations. volume 63, pages 608–616. Elsevier, 2012.

[5] Prakash Kumar Sahu and S Saha Ray. Chebyshev wavelet method for numeri-
cal solutions of integro-differential form of lane–emden type differential equations.



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

International Journal of Wavelets, Multiresolution and Information Processing,
15(02):1750015, 2017.

[6] Nidal Anakira, Ali Jameel, Mohmmad Hijazi, Abedel-Karrem Alomari, and Noraziah
Man. A new approach for solving multi-pantograph type delay differential equations.
International Journal of Electrical & Computer Engineering (2088-8708), 12(2), 2022.

[7] Abdullah Dawar, Hamid Khan, Saeed Islam, and Waris Khan. The improved residual
power series method for a system of differential equations: a new semi-numerical
method. International Journal of Modelling and Simulation, 45(4):1272–1285, 2025.

[8] Dulfikar Jawad Hashim, Ali Fareed Jameel, Teh Yuan Ying, AK Alomari, and
NR Anakira. Optimal homotopy asymptotic method for solving several models of
first order fuzzy fractional ivps. volume 61, pages 4931–4943. Elsevier, 2022.

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

[10] 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.

[11] Rashid Nawaz, Abraiz Khattak, Muhammad Akbar, Sumbal Ahsan, Zahir Shah, and
Adam Khan. Solution of fractional-order integro-differential equations using optimal
homotopy asymptotic method: R. nawaz et al. Journal of Thermal Analysis and
Calorimetry, 146(3):1421–1433, 2021.

[12] 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.

[13] M Matinfar, M Saeidy, and J Vahidi. Application of homotopy analysis method for
solving systems of volterra integral equations. Advances in Applied Mathematics and
Mechanics, 4(1):36–45, 2012.

[14] Ahmed A Hamoud, LAFTA A Dawood, KP Ghadle, and SM Atshan. Usage of
the modified variational iteration technique for solving fredholm integro-differential
equations. volume 9, pages 895–902, 2019.

[15] Zhimin Hong, Xiangzhong Fang, Zaizai Yan, and Hui Hao. On solving a system
of volterra integral equations with relaxed monte carlo method. Journal of Applied
Mathematics and Physics, 4(7):1315–1320, 2016.

[16] Praveen Agarwal, Abd-Allah Hyder, Mohammed Zakarya, Ghada AlNemer, Clemente
Cesarano, and Dario Assante. Exact solutions for a class of wick-type stochastic (3+
1)-dimensional modified benjamin–bona–mahony equations. Axioms, 8(4):134, 2019.

[17] 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.

[18] N Rabit Anakira, AK Alomari, AF Jameel, and Ishak Hashim. Multistage optimal
homotopy asymptotic method for solving initial-value problems. J. Nonlinear Sci.



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

Appl, 9(4):1826–1843, 2016.
[19] 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.

[20] Abdul Hadi Bhatti, Sharmila Karim, Ala Amourah, Ali Fareed Jameel, Feras Yousef,
and Nidal Anakira. A novel approximation method for solving ordinary differential
equations using the representation of ball curves. Mathematics, 13(2):250, 2025.

[21] Nidal Anakira, Osama Oqilat, Adel Almalki, Irianto Irianto, Saad Meqdad, and Ala
Amourah. Numerical proceduers for computing the exact solutions to systems of
ordinary differential equations. Int. J. Neutrosophic Sci, 2:165–65, 2025.

[22] AK Alomari, N Ratib Anakira, and I Hashim. Multiple solutions of problems in
fluid mechanics by predictor optimal homotopy asymptotic method. Advances in
Mechanical Engineering, 6:372537, 2014.

[23] Ala Amourah, Omar Alnajar, Maslina Darus, Ala Shdouh, and Osama Ogilat. Esti-
mates for the coefficients of subclasses defined by the bell distribution of bi-univalent
functions subordinate to gegenbauer polynomials. Mathematics, 11(8):1799, 2023.

[24] Abdullah Alsoboh, Ala Amourah, Maslina Darus, and Carla Amoi Rudder. Investi-
gating new subclasses of bi-univalent functions associated with q-pascal distribution
series using the subordination principle. Symmetry, 15(5):1109, 2023.

[25] 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.

[26] N Ratib Anakira, AK Alomari, and Ishak Hashim. Application of optimal homotopy
asymptotic method for solving linear delay differential equations. 1571(1):1013–1019,
2013.

[27] Sehar Iqbal, M Idrees, Abdul Majeed Siddiqui, and Ali R Ansari. Some solutions of the
linear and nonlinear klein–gordon equations using the optimal homotopy asymptotic
method. Applied mathematics and computation, 216(10):2898–2909, 2010.

[28] Mohsen Sheikholeslami, Hamid Reza Ashorynejad, Davood Domairry, Ishak Hashim,
et al. Investigation of the laminar viscous flow in a semi-porous channel in the presence
of uniform magnetic field using optimal homotopy asymptotic method. volume 41,
pages 1281–1285, 2012.

[29] Muhammad Sadiq Hashmi, Nasir Khan, and Sehar Iqbal. Optimal homotopy asymp-
totic method for solving nonlinear fredholm integral equations of second kind. Applied
Mathematics and Computation, 218(22):10982–10989, 2012.

[30] NR Anakira, AK Alomari, AF Jameel, and I Hashim. Multistage optimal homo-
topy asymptotic method for solving boundary value problems with robin boundary
conditions. Far East Journal of Mathematical Sciences, 102(8):1727–1744, 2017.

[31] Rashid Nawaz, Sumbal Ahsan, Muhammad Akbar, M Farooq, M Sulaiman, Hakeem
Ullah, and Saeed Islam. Semi analytical solutions of second type of three-dimensional
volterra integral equations. International Journal of Applied and Computational
Mathematics, 6(4):109, 2020.

[32] Rashid Nawaz, Abraiz Khattak, Muhammad Akbar, Sumbal Ahsan, Zahir Shah, and



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

Adam Khan. Solution of fractional-order integro-differential equations using optimal
homotopy asymptotic method: R. nawaz et al. Journal of Thermal Analysis and
Calorimetry, 146(3):1421–1433, 2021.