J. Nig. Soc. Phys. Sci. 3 (2021) 159–164

Journal of the
Nigerian Society

of Physical
Sciences

A New Multi-Step Method for Solving Delay Differential
Equations using Lagrange Interpolation

V. J. Shaalinia, S. E. Fadugbab,∗

aDepartment of Mathematics, Bishop Heber College, Trichy, India
bDepartment of Mathematics, Ekiti State University, Ado Ekiti, Nigeria

Abstract

This paper presents 2-step p-th order ( p = 2, 3, 4) multi-step methods that are based on the combination of both polynomial and exponential func-
tions for the solution of Delay Differential Equations (DDEs). Furthermore, the delay argument is approximated using the Lagrange interpolation.
The local truncation errors and stability polynomials for each order are derived. The Local Grid Search Algorithm (LGSA) is used to determine
the stability regions of the method. Moreover, applicability and suitability of the method have been demonstrated by some numerical examples of
DDEs with constant delay, time dependent and state dependent delays. The numerical results are compared with the theoretical solution as well
as the existing Rational Multi-step Method2 (RMM2).

DOI:10.46481/jnsps.2021.247

Keywords: Multi-step method, Delay differential equations, Interpolating function, Lagrange interpolation, Stability polynomial, Stability region

Article History :
Received: 08 June 2021
Received in revised form: 12 July 2021
Accepted for publication: 24 July 2021
Published: 29 August 2021

c©2021 Journal of the Nigerian Society of Physical Sciences. All rights reserved.
Communicated by: T. Latunde

1. Introduction

Delay Differential Equations (DDEs) are differential equa-
tions in which the derivative of the unknown function depends
not only at its present time but also at the previous times. In
Ordinary Differential Equations (ODEs), a simple initial con-
dition is given. But to specify DDEs, additional information
is needed. Because the derivative depends on the solution at
the previous times, an initial history function which gives in-
formation about the solution in the past needs to be specified.
A general form of the first order DDE is

y′ (t) = f (t, y (t) , y (t −τ)) , t > t0

∗Corresponding author tel. no: +2348067032044
Email address: sunday.fadugba@eksu.edu.ng (S. E. Fadugba )

y (t) = ϑ (t) , t ≤ t0 (1)

where ϑ(t) is the initial function and τ is the delay term. The
function ϑ(t) is also known as the ‘history function’, as it gives
information about the solution in the past. If the delay term τ
is a constant, then it is called constant delay. If it is function of
time t, then it is called time dependent delay. If it is a function
of time t and y(t), then it is called state dependent delay.

These equations arise in population dynamics, control sys-
tems, chemical kinetic, and in several areas of science and en-
gineering [1, 2, 3]. Recently there has been a growing interest
in obtaining the numerical solutions of DDEs. Rostann et al.
[4] implemented Adomian decomposition method for the solu-
tion of system of DDEs. Two and three point one-step block
method for solving DDEs was developed by [5]. Block method
for solving Pantograph type functional DDEs was described by

159



Shaalini & Fadugba / J. Nig. Soc. Phys. Sci. 3 (2021) 159–164 160

[6]. An exact/approximate solution of DDEs by using the com-
bination of Laplace and the variational iteration method were
obtained by [7]. RK method based on Harmonic Mean for solv-
ing DDEs with constant lags was proposed by [8]. Several nu-
merical methods have been constructed for solving stiff DDEs,
see [9, 10, 11]. Several multi-step techniques using varieties of
interpolating polynomials and functions have been developed
to solve ODEs such as [12, 13, 14, 15, 16, 17], just to mention
a few.

In this paper, we present the 2-step p-th order ( p = 2, 3, 4)
multi-step method for solving DDEs. This method has been
referred here as EPMM (2, p), (p = 2, 3, 4). The organization
of this paper is as follows: In Section Two, the derivation of
EPMM (2, p), (p = 2, 3, 4) is given. In Section Three, the sta-
bility analysis of EPMM (2, p), (p = 2, 3, 4) has been presented.
In Section Four, numerical illustrations of DDEs are provided.
Moreover, the numerical results are compared with the exist-
ing Rational Multi-step Method2 (RMM2) to demonstrate the
efficiency and suitability of the method.

2. Derivation of EPMM (2, p)

For 2-step p-th order EPMM, let us assume an approxima-
tion to the analytical solution y (tn+2) of (1) given by

yn+2 = a0e
2h

+ 1 +
p∑

j=1

b jh
j, (2)

where a0, b j, ( j=1,2,. . . p) are parameters that may contain the
approximation of y (tn) and higher derivatives of y (tn).
With EPMM in (2) we associate the difference operator L de-
fined by

L[y (t) ; h]EPM M =

y (t + 2h) − (1 + p∑
j=0

b jh
j)

 − a0e2h(3)
where y(t) is an arbitrary, continuous and differentiable func-
tion.
Expanding y(t + 2h) as Taylor series and collecting terms in (3),

L[y (t) ; h]EPM M = C0h
0
+C1h

1
+· · ·+C ph

p
+C p+1h

p+1
+. . .(4)

where Ci, i = 0, 1, . . . , p, p + 1 are the coefficients that need
to be determined.
For 2-step second order EPMM, we take p = 2 and expand y(t +
2h) via Taylor series, (3) becomes

L
[
y (t) ; h

]
EPM M(2,2) = −1 + y (t) + h

(
−b1 + 2y

′

(t)
)

+h2
(
2y
′′

(t) − b2
)

+ h3
(

4
3

y
′′′

(t)
)

+ a0e
2h

+ O(h4) (5)

Using e2h ≈ 1 + 2h + 2h2 in (2), we get

L
[
y (t) ; h

]
EPM M(2,2) = −1 − a0 + y (t) + h

(
−b1 − 2a0 + 2y

′

(t)
)

+h2
(
2y
′′

(t) − 2a0 − b2
)
+ h3

(
4
3

y
′′′

(t)
)

+ a0e
2h

+ O(h4)(6)

Comparing (5) and (6), we have
C0 = −(1 + a0) + y (t) ,
C1 = −b1 − 2a0 + 2y

′

(t) ,
C2 = 2y

′′

(t) − 2a0 − b2,
C3 =

4
3 y
′′′

(t)

(7)

For second order EPMM, we put C0 = C1 = C2 = 0 in (7) and
get the following solutions:

a0 = y (t)−1, b1 = 2(y
′

(t)−y(t)+1), b2 = 2(y
′′

(t)−y (t)+1)(8)

If we write yn = y(tn) and y
(m)
n = y

(m)(tn) for m = 1, 2,. . . , then
(8) becomes

A = yn,
b1 = 2

(
yn
′

− yn + 1
)
,

b2 = 2(yn
′′

− yn + 1)
(9)

Taking p = 2 and e2h ≈ 1 + 2h + 2h2 in (2), we get

yn+2 = (a0 + 1) + h (2a0 + b1) + h
2(2a0 + b2) (10)

Substituting (9) into (10), we have

yn+2 = yn + 2hyn
′

+ 2h2yn
′′

(11)

The local truncation error of EPMM (2, 2) is given by,

LT E EPM M(2,2) = h
3
(

4
3

yn
′′′

)
+ O(h4)

Taking p = 3 in (2) and on simplification, we get the formula
for EPMM (2, 3)

yn+2 = yn + 2hyn
′

+ 2h2yn
′′

+
4
3

h3yn
′′′

(12)

The local truncation error of EPMM (2, 3) is given by,

LT E EPM M(2,3) = h
4
(

2
3

yn
(4)

)
+ O(h5)

Taking p = 4 in (2) and on simplification, we get the formula
for EPMM (2, 4)

yn+2 = y (t)+2hy
′

(t)+2h2y
′′

(t)+
4
3

h3y
′′′

(t)+
2
3

h4y(4) (t)(13)

The local truncation error of EPMM (2, 4) is given by,

LT E EPM M(2,4) = h
5
(

4
15

yn
(5)

)
+ O(h6)

3. Stability Analysis of EPMM

In this section, we derive the stability polynomials of EPMM
(2, p), (p = 2, 3, 4) and their corresponding stability regions
were obtained.
We consider a commonly used linear test equation with a con-
stant delay τ = mh where m is a positive integer,

y
′

(t) = λy (t) + µy(t −τ), t > t0
160



Shaalini & Fadugba / J. Nig. Soc. Phys. Sci. 3 (2021) 159–164 161

y (t) = φ(t), t ≤ t0 (14)

where λ,µ ∈ C, τ > 0 andΦ is continuous.
Using (11) in (14), we get

yn+2 = yn+2h (λyn + µy (tn −τ))+2h
2 (λy′n + µy′ (tn −τ))(15)

y (tn − mh) = y (tn−m) =
s1∑

l=−r1

Ll (ci)yn−m+l

with

Ll (ci) =
s1∏

j=−r1

ci − j1
l − j1

, j1 , l and r1, s1 > 0

Taking

y (tn −τ) =
s1∑

l=−r1

Ll(c)yn−m+l

and

y
′

(tn −τ) = λ
s1∑

l=−r1

Ll (c) yn−m+l + µ
s1∑

l=−r1

Ll (c) yn−2m+l (16)

Then (15) becomes

yn+2 = yn + 2h

λyn + µ s1∑
l=−r1

Ll (c) yn−m+l


+2h2


λ
(
λyn + µ

∑s1
l=−r1

Ll (c) yn−m+l
)

+µλ
∑s1

l=−r1
Ll (c) yn−m+l

+µ
∑s1

l=−r1
Ll (c) yn−2m+l


yn+2 = yn + 2λhyn + 2λ

2h
2
yn

+

s1∑
l=−r1

Ll (c) yn−m+l
(
2µh + 4h2µλ

)
+

s1∑
l=−r1

Ll (c) yn−2m+l
(
2h2µ2

)

yn+2 = yn
(
1 + 2λh + 2(λh)2

)
+

s1∑
l=−r1

Ll (c) yn−m+l (µh (2 + 4λh))

+

s1∑
l=−r1

Ll (c) yn−2m+l
(
2(µh)2

)
Let α = λh and β = µh then the above equation becomes

yn+2 = yn
(
1 + 2α + 2α2

)
+ (β (2 + 4α))

s1∑
l=−r1

Ll (c) yn−m+l

+2β2
s1∑

l=−r1

Ll (c) yn−2m+l

To obtain the stability polynomial, the delay term is approxi-
mated using three points Lagrange interpolation.
By putting n − m + l = 0 and n − 2m + l = 0 and by taking l =
-1, 0, 1, the stability polynomial will be in the standard form.

The recurrence is stable if the zeros of ζi of the stability poly-
nomial

S (α,β : ζ) = ζn+2 −
(
1 + α + 2α2

)
ζn

−β (2 + 4α)
(
L−1 (c) + L0 (c) ζ + L1 (c) ζ

2
)

−2β2
(
L−1 (c) + L0 (c) ζ + L1 (c) ζ

2
)

satisfies the root condition|ζi| ≤ 1. From this, the stability poly-
nomial for the method GRMM (2, 2) with τ = 1 is given as

S (α,β : ζ) = ζn+2 −
(
1 + α + 2α2

)
ζn −

(
2β + 2β2 + 4αβ

)
Similarly, by considering suitable number of points in Lagrange
interpolation according to the order of the method, we can ob-
tain the corresponding stability polynomials of EPMM (2, p).
When p = 3, the stability polynomial for EPMM (2, 3) is given
as

S (α,β : ζ) = ζn+2 −
(
1 + α + 2α2 +

4
3
α3

)
ζn

−

(
2β + 2β2 +

4
3
β3 + 4αβ + 4α2β + 4αβ2

)
When p = 4, the stability polynomial for EPMM (2, 4) is given
as

S (α,β : ζ) = ζn+2 −
(
1 + α + 2α2 +

4
3
α3 +

2
3
α4

)
ζn

−

(
2β + 2β2 +

4
3
β3 +

2
3
β4 + 2αβ + 4α2β + 4αβ2 +

8
3
αβ3

+
8
3
α3β + 4α2β2

)
The stability regions of EPMM (2, 2), EPMM (2, 3) and

EPMM (2, 4) are given in Figures 1 -3.
In a similar manner, we can obtain the stability polynomials
and their corresponding regions of EPMM with r-step and of
any order p.

4. Numerical Examples

Example 1: (Stiff linear system with multiple delays)

y1
′

(t) = −
1
2

y1 (t) −
1
2

y2 (t − 1) + f1 (t) ,

y2
′

(t) = −y2 (t) −
1
2

y1

(
t −

1
2

)
+ f2(t), 0 ≤ t ≤ 1

with initial conditions

y1 (t) = e
−t/2,

−1
2
≤ t ≤ 0,

y2 (t) = e
−t, −1 ≤ t ≤ 0

161



Shaalini & Fadugba / J. Nig. Soc. Phys. Sci. 3 (2021) 159–164 162

Table 1. Comparison of Absolute Error Results in EPMM and RMM2 for Example 1
Time (t) Y EPMM (2, 2) RMM2 (2, 2) EPMM (2, 3) RMM2 (2,3) EPMM (2, 4) RMM2(2, 4)

0.2 y1 1.52e-06 7.54e-07 3.80e-09 1.26e-09 7.60e-12 9.06e-07
y2 3.26e-04 3.31e-04 3.26e-04 3.26e-04 3.15e-04 3.32e-04

0.4 y1 2.75e-06 1.36e-06 6.88e-09 2.28e-09 1.38e-11 1.64e-06
y2 5.62e-04 5.71e-04 5.62e-04 5.62e-04 5.43e-04 5.72e-04

0.6 y1 3.73e-06 1.85e-06 9.33e-09 3.10e-09 1.87e-11 2.22e-06
y2 7.27e-04 7.38e-04 7.27e-04 7.27e-04 7.04e-04 7.40e-04

0.8 y1 4.50e-06 2.23e-06 1.13e-08 3.73e-09 2.25e-11 2.68e-06
y2 8.36e-04 8.48e-04 8.36e-04 8.36e-04 8.12e-04 8.51e-04

1.0 y1 5.09e-06 2.53e-06 1.27e-08 4.22e-09 2.55e-11 3.04e-06
y2 9.03e-04 9.16e-04 9.03e-04 9.03e-04 8.78e-04 9.18e-04

Table 2. Comparison of Absolute Error Results in EPMM and RMM2 for Example 2
Time (t) EPMM (2, 2) RMM2 (2, 2) EPMM (2, 3) RMM2 (2, 3) EPMM (2, 4) RMM2 (2, 4)

1.1 1.82e-06 4.04e-06 1.31e-06 1.46e-06 1.35e-06 9.58e-06
1.2 9.26e-07 2.51e-06 1.25e-06 1.20e-06 2.75e-06 8.03e-06
1.3 1.50e-06 3.23e-08 2.22e-06 4.13e-06 2.26e-07 8.41e-06
1.4 3.50e-06 5.92e-06 2.70e-06 3.29e-06 2.99e-06 1.19e-05
1.5 4.95e-05 1.96e-06 3.13e-06 6.16e-06 4.39e-06 6.83e-06

Table 3. Comparison of Absolute Error Results in EPMM and RMM2 for Example 3
Time(t) EPMM (2, 2) RMM2 (2, 2) EPMM (2, 3) RMM2 (2, 3) EPMM (2, 4) RMM2 (2, 4)

0.2 1.33e-05 1.35e-05 1.97e-09 2.36e-07 4.92e-09 4.87e-07
0.4 2.60e-05 2.80e-05 2.46e-08 5.62e-07 1.76e-09 5.46e-06
0.6 3.78e-05 4.48e-05 5.67e-08 7.45e-07 3.03e-09 3.58e-05
0.8 4.79e-05 6.56e-05 9.93e-08 8.72e-07 2.45e-09 2.32e-04
1.0 5.63e-05 9.31e-05 2.85e-06 2.85e-06 3.24e-09 3.24e-04

Figure 1. Stability Region of 2-step Second Order EPMM

where
f1 (t) =

1
2

e−(t−1)

and
f2 (t) =

1
2

e−(t−1/2)/2

The exact solution is given by

y1 (t) = e
−t/2, y2 (t) = e

−t

Figure 2. Stability Region of 2-step Third Order EPMM

Example 2: (Time-dependent delay)

y
′

(t) =
t − 1

t
y(ln (t) − 1)y(t), 1 ≤ t ≤

3
2

With initial condition

y (t) = 1, 0 ≤ t ≤ 1

162



Shaalini & Fadugba / J. Nig. Soc. Phys. Sci. 3 (2021) 159–164 163

Figure 3. Stability Region of 2-step Fourth Order EPMM

Figure 4. Comparison of Error Graph of y1 and y2 in Example 1

Figure 5. Comparison of Absolute Error Graph of y in Example 2

and the exact solution is given by

y (t) = exp(t − ln (t) − 1), 1 ≤ t ≤
3
2

Example 3: (State-dependent delay)

y
′

(t) = cos(t)y
(
y (t)−2

)
, t ≥ 0

Figure 6. Comparison of Absolute Error Graph of y in Example 3

With initial condition,

y (t) = 1, t ≤ 0

and the exact solution is given by

y (t) = sin (t) + 1, 0 ≤ t ≤ 1

By taking the step-size h = 0.01in the above examples, the
absolute errors by using EPMM and RMM2 are given in Tables
1 – 3 and their corresponding error graphs are shown in Figures
4 – 6.

5. Conclusion

In this paper, the new multi-step method of r-step and p-th
order that are based on interpolating functions which consists of
both polynomial and exponential function is presented for solv-
ing DDEs. The local truncation errors have been determined.
The stability polynomials of EPMM (2, p) where p = 2, 3, 4 are
derived and their corresponding stability regions are obtained
and shown in Figures 1–3. The delay argument is approximated
using Lagrange interpolation. Numerical examples of DDEs
with constant delay, time dependent delay and state dependent
delays have been considered to demonstrate the efficiency of
the proposed method. The comparative absolute error analyses
of EPMM (2,p) in the context of RMM2 (2,p) for Examples 1,
2 and 3 were shown in Tables 1, 2 and 3, respectively. From the
Figures 4 – 6, it is evident that the newly proposed method gives
results with good accuracy than the existing RMM2. Hence, it
is concluded that the proposed EPMM (2, p) is suitable for solv-
ing DDEs.

Acknowledgments

We thank the referees and the editor for their contributions
on the improvement of this paper.

163



Shaalini & Fadugba / J. Nig. Soc. Phys. Sci. 3 (2021) 159–164 164

References

[1] Y. Kuang, Delay differential equations with applications in population
dynamics, Academic Press, Boston, San Diego, New York. (1993).

[2] E. Fridman, L. Fridman & E. Shusti, “Steady modes in relay control sys-
tems with time delay and periodic disturbances”, Journal of Dynamical
Systems Measurement and Control 122 (2000) 732.

[3] I. Epstein & Y. Luo, “Differential delay equations in chemical kinetics:
non-linear models; the cross-shaped phase diagram and the oregonator”,
Journal of Chemical Physic 95 (1991) 244.

[4] R. K. Saeed & B. M. Rahman, “Adomian decomposition method for solv-
ing system of delay differential equation”, Australian Journal of Basic and
Applied Sciences 4 (2010) 3613.

[5] H. M. Radzi, Z. A. Majid, F. Ismail & M. Suleiman, “Two and three point
one-step block method for solving delay differential equations”, Journal
of Quality Measurement and Analysis 8 (2012) 29.

[6] F. Ishak, M. B. Suleiman & Zanariah A. Majid, “Block Method for solv-
ing Pantograph-type functional differential equations”, Proceeding of the
World Congress on Engineering 2 (2013).

[7] T. A. Biala, O. O. Asim & Y. O. Afolabi, “A combination of the Laplace
and the variational iteration method for the analytical treatment of delay
differential equations”, International Journal of Differential Equations and
Applications 13 (2014) 164.

[8] A. E. K. Puhpam & J. V. Shaalini, “Solving delay differential equations
with constant lags using RKHaM method”, International Journal of Sci-
entific Research 5 (2016) 585.

[9] G. A. Bocharov, G. I. Marchuk & A .A. Romanyukha, “Numerical so-

lution by LMMs of stiff delay differential systems modelling an immune
response”, Numerische Mathematik 73 (1996) 131.

[10] J. Vinci Shaalini & A. Emimal Kanaga Pushpam, “Analysis of Composite
Runge Kutta Methods and New One-step Technique for stiff Delay Differ-
ential Equations”, IAENG International Journal of Applied Mathematics
4 (2019) 359.

[11] J. V. Shaalini & A. E. K. Pushpam, “A new one step method for solving
stiff and non-stiff delay differential equations using Lagrange interpola-
tion”, Journal of Applied Science and Computations 6 (2019) 949.

[12] S. O. Fatunla, “Nonlinear multistep methods for initial value problems”,
Computers and Mathematics with Applications 8 (1982) 231.

[13] K. O. Okosun & R. A. Ademiluyi, “A two-step second order inverse poly-
nomial methods for integration of differential equations with singulari-
ties”, Research Journal of Applied Sciences 2 (2007) 13.

[14] O. E. Abolarin & S. W. Akingbade, “Derivation and application of fourth
stage inverse polynomial scheme to initial value problems”, IAENG In-
ternational Journal of Applied Mathematics 47 (2017) 459.

[15] T. Y. Ying & N. Yaacob, “A new class of rational multi-step methods
for solving initial value problems”, Malaysian Journal of Mathematical
Sciences, 7 (2013) 31.

[16] J. Vinci Shaalini, A. & E. K. Pushpam, “A Class of 2-Step fourth order ra-
tional multi-step methods for solving delay differential equations”, Inter-
national Journal for Research in Engineering Application & Management
5 (2019) 164.

[17] S. E. Fadugba, S.N. Ogunyebi & B.O. Falodun, “An examination of a sec-
ond order numerical method for solving initial value problems”, Journal
of the Nigerian Society of Physical Sciences 2 (2020) 120.

164