Physics Letters A 372 (2008) 470481

www.elsevier.com/locate/pla
Adaptation of homotopy-perturbation method for numericanalytic
solution of system of ODEs
I. Hashim

, M.S.H. Chowdhury
School of Mathematical Sciences, Universiti Kebangsaan Malaysia, 43600 Bangi, Selangor, Malaysia
Received 21 April 2007; received in revised form 11 July 2007; accepted 26 July 2007
Available online 2 August 2007
Communicated by A.R. Bishop
Abstract
A new reliable algorithm based on an adaptation of the standard Homotopy-Perturbation Method (HPM) is presented. The HPM is treated,
for the rst time, as an algorithm in a sequence of intervals (i.e., time step) for nding accurate approximate solutions of linear and nonlinear
systems of ODEs. Numerical comparisons between the Multistage Homotopy-Perturbation Method (MHPM) and the available exact solution and
the classical fourth-order RungeKutta (RK4) method reveal that the new technique is a promising tool for linear and nonlinear systems of ODEs.
2007 Elsevier B.V. All rights reserved.
Keywords: Homotopy-perturbation method; RungeKutta method; System of ODEs
1. Introduction
It is well known that many physical and engineering phenomena can be modelled by systems of ODEs. One important class of
these systems is the general system of rst-order ODEs:
(1)
du
1
dt
+g
1
(t, u
1
, u
2
, . . . , u
m
) =f
1
(t ),
(2)
du
2
dt
+g
2
(t, u
1
, u
2
, . . . , u
m
) =f
2
(t ),
.
.
.
(3)
du
m
dt
+g
m
(t, u
1
, u
2
, . . . , u
m
) =f
m
(t ),
subject to the initial conditions
(4) u
1
(t
0
) = c
1
, u
2
(t
0
) =c
2
, . . . , u
m
(t
0
) =c
m
.
Finding accurate and efcient methods for solving (1)(4) has long been an active research undertaking. Some classes of system
(1)(4) have been solved by the numericanalytic technique of the Adomian Decomposition Method (ADM) in [17]. However,
one notable difculty inherent in ADM is the calculation of the so-called Adomian polynomials which can be cumbersome in
general. Another analytic method which has been shown to be much simpler than the ADM is called the homotopy-perturbation
method (HPM), rst developed by He [815]. The method yields rapidly convergent series solutions [1620]. Very recently, HPM
*
Corresponding author.
E-mail address: ishak_h@ukm.my (I. Hashim).
0375-9601/$ see front matter 2007 Elsevier B.V. All rights reserved.
doi:10.1016/j.physleta.2007.07.067
I. Hashim, M.S.H. Chowdhury / Physics Letters A 372 (2008) 470481 471
was employed for solving singular second-order differential equations [21], nonlinear population dynamics models [22] and time-
dependent EmdenFowler type equations [23].
Like ADM and based on our observations, it is hopeless to nd HPM solutions of system (1)(4) valid globally in time. This
lack of global convergence is overcome, in this Letter, by recursively applying HPM over successive time intervals, thereby solving
a sequence of IVPs, each valid over a given time interval. We shall call this numericanalytic adaptation of HPM, the multistage
HPM (MHPM). Numerical comparisons with the exact, HPM and fourth-order RungeKutta (RK4) solutions show that MHPM is
accurate.
2. HPM solutions
First write system (1)(3) in the operator form:
(5) L(u
1
) +N
1
(u
1
, u
2
, . . . , u
m
) f
1
=0,
(6) L(u
2
) +N
2
(u
1
, u
2
, . . . , u
m
) f
2
=0,
.
.
.
(7) L(u
m
) +N
m
(u
1
, u
2
, . . . , u
m
) f
m
=0,
subject to the initial conditions (4), where L = d/dt is linear operator and N
1
, N
2
, . . . , N
m
are nonlinear operators. We shall next
present the solution approaches of (5)(7) based on the standard HPM and MHPM separately.
2.1. Solution by standard HPM
According to HPM, we construct a homotopy for (5)(7) which satises the following relations:
(8) L(u
1
) L(v
1
) +pL(v
1
) +p
_
N
1
(u
1
, u
2
, . . . , u
m
) f
1
_
=0,
(9) L(u
2
) L(v
2
) +pL(v
2
) +p
_
N
2
(u
1
, u
2
, . . . , u
m
) f
2
_
=0,
.
.
.
(10) L(u
m
) L(v
m
) +pL(v
m
) +p
_
N
m
(u
1
, u
2
, . . . , u
m
) f
m
_
=0,
where p [0, 1] is an embedding parameter and v
1
, v
2
, . . . , v
m
are initial approximations which satisfying the given conditions. It
is obvious that when the perturbation parameter p = 0, Eqs. (8)(10) become a linear system of equations and when p = 1 we get
the original nonlinear system of equations.
Let us take the initial approximations as follows:
(11) u
1,0
(t ) =v
1
(t ) =u
1
(t
0
) =c
1
,
(12) u
2,0
(t ) =v
2
(t ) =u
2
(t
0
) =c
2
,
.
.
.
(13) u
m,0
(t ) =v
m
(t ) =u
m
(t
0
) =c
m
,
and
(14) u
1
(t ) =u
1,0
(t ) +pu
1,1
(t ) +p
2
u
1,2
(t ) +p
3
u
1,3
(t ) + ,
(15) u
2
(t ) =u
2,0
(t ) +pu
2,1
(t ) +p
2
u
2,2
(t ) +p
3
u
2,3
(t ) + ,
.
.
.
(16) u
m
(t ) =u
m,0
(t ) +pu
m,1
(t ) +p
2
u
m,2
(t ) +p
3
u
m,3
(t ) + ,
where u
i,j
(i =1, 2, . . . , m; j =1, 2, . . .) are functions yet to be determined. Substituting (11)(16) into (8)(10) and arranging the
coefcients of the same powers of p, we get
(17) L(u
1,1
) +L(v
1
) +N
1
(u
1,0
, u
2,0
, . . . , u
m,0
) f
1
=0, u
1,1
(t
0
) =0,
(18) L(u
2,1
) +L(v
2
) +N
2
(u
1,0
, u
2,0
, . . . , u
m,0
) f
2
=0, u
2,1
(t
0
) =0,
.
.
.
(19) L(u
m,1
) +L(v
m
) +N
m
(u
1,0
, u
2,0
, . . . , u
m,0
) f
m
=0, u
m,1
(t
0
) =0,
472 I. Hashim, M.S.H. Chowdhury / Physics Letters A 372 (2008) 470481
(20) L(u
1,2
) +N
1
(u
1,1
, u
2,1
, . . . , u
m,1
) =0, u
1,2
(t
0
) =0,
(21) L(u
2,2
) +N
2
(u
1,1
, u
2,1
, . . . , u
m,1
) =0, u
2,2
(t
0
) =0,
.
.
.
(22) L(u
m,2
) +N
m
(u
1,1
, u
2,1
, . . . , u
m,1
) =0, u
m,2
(t
0
) =0,
etc. We solve the above systems of equations for the unknowns u
i,j
(i =1, 2, . . . , m; j =1, 2, . . .) by applying the inverse operator
(23) L
1
() =
t
_
0
() dt.
Therefore, according to HPM the n-term approximations for the solutions of (5)(7) can be expressed as
(24)
1,n
(t ) =u
1
(t ) = lim
p1
u
1
(t ) =
n1

k=0
u
1,k
(t ),
(25)
2,n
(t ) =u
2
(t ) = lim
p1
u
2
(t ) =
n1

k=0
u
2,k
(t ),
.
.
.
(26)
m,n
(t ) =u
m
(t ) = lim
p1
u
m
(t ) =
n1

k=0
u
m,k
(t ).
2.2. Solution by MHPM
The approximate solutions (24)(26) are generally, as shall be shown in the numerical experiments of this Letter, not valid for
large t . A simple way of ensuring validity of the approximations for large t is to treat (17)(22) as an algorithm for approximating
the solutions of (1)(4) in a sequence of intervals choosing the initial approximations as
(27) u
1,0
(t ) =v
1
(t ) =u
1
(t

) =c

1
,
(28) u
2,0
(t ) =v
2
(t ) =u
2
(t

) =c

2
,
.
.
.
(29) u
m,0
(t ) =v
m
(t ) =u
m
(t

) =c

m
.
Now we solve (17)(22) for the unknowns u
i,j
(i =1, 2, . . . , m; j =1, 2, . . .) by applying the inverse linear operator
(30) L
1
() =
t
_
t

() dt.
In order to carry out the iterations in every subinterval of equal length t , [0, t
1
), [t
1
, t
2
), [t
2
, t
3
), . . . , [t
j1
, t ), we would need to
know the values of the following,
(31) u

1,0
(t ) =u
1
(t

), u

2,0
(t ) =u
2
(t

), . . . , u

m,0
(t ) =u
m
(t

).
But, in general, we do not have these information at our clearance except at the initial point t

= t
0
. A simple way for obtaining the
necessary values could be by means of the previous n-term approximations
1,n
,
2,n
, . . . ,
m,n
of the preceding subinterval, i.e.,
(32) u

1,0

1,n
(t

), u

2,0

2,n
(t

), . . . , u

m,0

m,n
(t

).
3. The fourth-order RungeKutta method (RK4)
In numerical analysis, the RungeKutta methods are an important family of implicit and explicit iterative methods for the
approximation of solutions of differential equations. One member of the family of RungeKutta methods is so commonly used is
the fourth-order RungeKutta method which is often referred as RK4. To illustrate the basic idea of RK4 method, consider an
initial value problem be specied as follows
(33) y

=f (t, y), y(t

0
) =y
0
.
I. Hashim, M.S.H. Chowdhury / Physics Letters A 372 (2008) 470481 473
Then, the RK4 method for this problem is given by the following equation:
(34) y
n+1
=y
n
+
h
6
(k
1
+2k
2
+2k
3
+k
4
),
where
(35) k
1
=f (t
n
, y
n
),
(36) k
2
=f
_
t
n
+
h
2
, y
n
+
h
2
k
1
_
,
(37) k
3
=f
_
t
n
+
h
2
, y
n
+
h
2
k
2
_
,
(38) k
4
=f (t
n
+h, y
n
+hk
3
).
Thus, the next value y
n+1
is determined by the present value y
n
plus the product of the size of the interval h and an estimated slope.
The slope is a weighted average of slopes:
k
1
is the slope at the beginning of the interval;
k
2
is the slope at the midpoint of the interval, using slope k
1
to determine the value of y at the point t
n
+ h/2 using Eulers
method;
k
3
is again the slope at the midpoint, but now using the slope k
2
to determine the y-value;
k
4
is the slope at the end of the interval, with its y-value determined using k
3
.
In averaging the four slopes, greater weight is given to the slopes at the midpoint:
(39) slope =
k
1
+2k
2
+2k
3
+k
4
6
.
The RK4 method has the error per step is on the order of h
5
, while the total accumulated error is of order h
4
.
4. Numerical experiments
In this section, we shall demonstrate the applicability and accuracy of MHPM to systems of linear and nonlinear ODEs through
three examples. The MHPM algorithm is coded in the computer algebra package Maple. We select a xed step sizes of t = 0.01
for the linear example and t = 0.001 for the nonlinear examples. The Maple environment variable digits controlling the number
of signicant digits is set to 16 in all the calculations done in this Letter.
4.1. Example 1
The rst system we shall study in this Letter is based on the simple linear system of ODEs considered in [24],
(40) x

=x +y,
(41) y

=x +y,
subject to the initial conditions
(42) x(0) =c
1
, y(0) =c
2
,
where the prime denotes differentiation with respect to t .
4.1.1. HPM solution
According to the HPM, we can construct a homotopy of system (40)(41) as follows:
(43) v

1
x

0
+p(x

0
v
1
v
2
) =0,
(44) v

2
y

0
+p(y

0
+v
1
v
2
) =0.
Let us choose the initial approximations as
(45) v
1,0
(t ) =x
0
(t ) =v
1
(0) =c
1
,
(46) v
2,0
(t ) =y
0
(t ) =v
2
(0) =c
2
,
474 I. Hashim, M.S.H. Chowdhury / Physics Letters A 372 (2008) 470481
and
(47) v
1
(t ) =v
1,0
(t ) +pv
1,1
(t ) +p
2
v
1,2
(t ) +p
3
v
1,3
(t ) + ,
(48) v
2
(t ) =v
2,0
(t ) +pv
2,1
(t ) +p
2
v
2,2
(t ) +p
3
v
2,3
(t ) + ,
where v
i,j
(i = 1, 2; j = 1, 2, 3, . . .) are functions yet to be determined. Substituting (45)(48) into (43)(44) and collecting terms
of the same powers of p, we have
(49) v

1,1
+x

0
v
1,0
v
2,0
=0, v
1,1
(0) =0,
(50) v

2,1
+y

0
+v
1,0
v
2,0
=0, v
2,1
(0) =0,
(51) v

1,2
v
1,1
v
2,1
=0, v
1,2
(0) =0,
(52) v

2,2
+v
1,1
v
2,1
=0, v
2,2
(0) =0,
(53) v

1,3
v
1,2
v
2,2
=0, v
1,3
(0) =0,
(54) v

2,3
+v
1,2
v
2,2
=0, v
2,3
(0) =0.
(55) v

1,4
v
1,3
v
2,3
=0, v
1,4
(0) =0,
(56) v

2,4
+v
1,3
v
2,3
=0, v
2,4
(0) =0,
(57) v

1,5
v
1,4
v
2,4
=0, v
1,5
(0) =0,
(58) v

2,5
+v
1,4
v
2,4
=0, v
2,5
(0) =0.
Solving the differential equations (49)(54) we get,
(59) v
1,1
(t ) =
t
_
0
v
1,0
dt +
t
_
0
v
2,0
dt =t c
1
+t c
2
,
(60) v
2,1
(t ) =
t
_
0
v
1,0
dt +
t
_
0
v
2,0
dt =c
1
t +c
2
t,
(61) v
1,2
(t ) =
t
_
0
v
1,1
dt +
t
_
0
v
2,1
dt =c
2
t
2
,
(62) v
2,2
(t ) =
t
_
0
v
1,1
dt +
t
_
0
v
2,1
dt =c
1
t
2
,
(63) v
1,3
(t ) =
t
_
0
v
1,2
dt +
t
_
0
v
2,2
dt =
1
3
c
2
t
3

1
3
c
1
t
3
,
(64) v
2,3
(t ) =
t
_
0
v
1,2
dt +
t
_
0
v
2,2
dt =
1
3
c
2
t
3

1
3
c
1
t
3
,
(65) v
1,4
(t ) =
t
_
0
v
1,3
dt +
t
_
0
v
2,3
dt =
1
6
c
1
t
4
,
(66) v
2,4
(t ) =
t
_
0
v
1,3
dt +
t
_
0
v
2,3
dt =
1
6
c
2
t
4
,
(67) v
1,5
(t ) =
t
_
0
v
1,4
dt +
t
_
0
v
2,4
dt =
1
30
c
1
t
5

1
30
c
2
t
5
,
(68) v
2,5
(t ) =
t
_
0
v
1,4
dt +
t
_
0
v
2,4
dt =
1
30
c
1
t
5

1
30
c
2
t
5
.
I. Hashim, M.S.H. Chowdhury / Physics Letters A 372 (2008) 470481 475
If 5-term approximations are sufcient, we obtain
(69) x(t ) = lim
p1
v
1
(t ) =
5

k=0
v
1,k
(t ),
(70) y(t ) = lim
p1
v
2
(t ) =
5

k=0
v
2,k
(t ),
which yield
(71) x(t ) =c
1
+c
1
t +c
2
t +c
2
t
2
+
1
3
c
2
t
3

1
3
c
1
t
3

1
6
c
1
t
4

1
30
c
1
t
5

1
30
c
2
t
5
,
(72) y(t ) =c
2
c
1
t +c
2
t c
1
t
2

1
3
c
2
t
3

1
3
c
1
t
3

1
6
c
2
t
4
+
1
30
c
1
t
5

1
30
c
2
t
5
.
Taking the actual physiological data [24,25], c
1
=0, c
2
=1, yields the 5-term HPM solutions,
(73) x(t ) =t +t
2
+
1
3
t
3

1
30
t
5
,
(74) y(t ) =1 +t
1
3
t
3

1
6
t
4

1
30
t
5
.
If more terms are added we approach the exact solution,
x(t ) =e
t
sint, y(t ) =e
t
cos t.
4.1.2. MHPM solution
Carrying out the steps involved in MHPM gives for all t t

,
(75) v
1,1
(t ) =
t
_
t

v
1,0
dt +
t
_
t

v
2,0
dt =(c

1
+c

2
)(t t

),
(76) v
2,1
(t ) =
t
_
t

v
1,0
dt +
t
_
t

v
2,0
dt =(c

1
+c

2
)(t t

),
(77) v
1,2
(t ) =
t
_
t

v
1,1
dt +
t
_
t

v
2,1
dt =c

2
(t t

)
2
,
(78) v
2,2
(t ) =
t
_
t

v
1,1
dt +
t
_
t

v
2,1
dt =c

1
(t t

)
2
,
(79) v
1,3
(t ) =
t
_
t

v
1,2
dt +
t
_
t

v
2,2
dt =
1
3
(c

2
c

1
)(t t

)
3
,
(80) v
2,3
(t ) =
t
_
t

v
1,2
dt +
t
_
t

v
2,2
dt =
1
3
(c

2
+c

1
)(t t

)
3
,
(81) v
1,4
(t ) =
t
_
t

v
1,3
dt +
t
_
t

v
2,3
dt =
1
6
c

1
(t t

)
4
,
(82) v
2,4
(t ) =
t
_
t

v
1,3
dt +
t
_
t

v
2,3
dt =
1
6
c

2
(t t

)
4
,
(83) v
1,5
(t ) =
t
_
t

v
1,4
dt +
t
_
t

v
2,4
dt =
1
30
(c

1
+c

2
)(t t

)
5
,
476 I. Hashim, M.S.H. Chowdhury / Physics Letters A 372 (2008) 470481
Table 1
Differences between the 5-term HPM and 5-term MHPM (t = 0.001) solutions with the exact solutions for the case c
1
=0 and c
2
= 1 related to example 1
=|Exact HPM| = |Exact MHPM|
t x y x y
0.05 1.749E-10 1.257E-12 5.780E-14 2.000E-15
0.5 1.859E-04 1.404E-05 7.976E-13 4.330E-13
1 1.264E-02 2.027E-03 1.627E-12 2.503E-12
1.5 1.514E-01 3.890E-02 5.550E-13 7.390E-12
2 8.816E-01 0.325E-01 6.747E-12 1.481E-11
2.5 0.341E+01 0.171E+01 2.720E-11 1.974E-11
3 1.007E+01 0.672E+01 6.820E-11 9.610E-12
3.5 2.415E+01 2.129E+01 1.102E-10 5.131E-11
4 4.852E+01 5.745E+01 1.718E-10 1.863E-10
4.5 8.161E+01 1.358E+02 1.324E-10 4.440E-10
5 1.098E+02 2.861E+02 1.232E-10 3.616E-10
(84) v
2,5
(t ) =
t
_
t

v
1,4
dt +
t
_
t

v
2,4
dt =
1
30
(c

1
c

2
)(t t

)
5
.
To carry out the iterations in every subinterval of equal length t , [0, t
1
), [t
1
, t
2
), [t
2
, t
3
), . . . , [t
n1
, t ), we would need to know the
values of the following,
(85) c

1
=x(t

), c

2
=y(t

).
In general, we do not have these information at our clearance except at the initial point t

= t
0
= 0 but we can obtain these values
following the MHPM as given in Section 2.2. We note that the 5-term approximations of x and y are denoted as x(t )
5
(t ) =

5
i=0
v
1,i
and y(t )
5
(t ) =

5
i=0
v
2,i
.
In Table 1 we present the differences between the 5-term HPM and 5-term MHPM solutions with the exact solutions for the case
c
1
=0 and c
2
=1. We observe that the 5-term MHPM solutions and the exact solutions agree with each other at least to 10-decimal
places, but the 5-term HPM solutions are valid only for very small time span t .
4.2. Example 2
The second system we shall study is a mathematical model of the beating action of the human heart given by the nonlinear
system of ODEs [24,25],
(86) x

=(T x y) x
3
,
(87) y

=x,
subject to the initial conditions
(88) x(0) =c
1
, y(0) =c
2
,
where x represents the length of muscle ber and y represents the amount of chemical control agent, and , T , c
1
are real constants.
4.2.1. HPM solution
According to the HPM, we can construct a homotopy of system (86)(87) as follows:
(89) v

1
x

0
+p
_
x

0
T v
1
+v
2
+v
3
1
_
=0,
(90) v

2
y

0
+p(y

0
v
1
) =0.
Let us choose the initial approximations as
(91) v
1,0
(t ) =x
0
(t ) =x(0) =c
1
,
(92) v
2,0
(t ) =y
0
(t ) =y(0) =c
2
.
Substituting (47)(48) and (91)(92) into (89) and (90), and collecting terms of the same powers of p, we have
(93) v

1,1
+x

0
T v
1,0
+v
2,0
+v
3
1,0
=0, v
1,1
(0) =0,
(94) v

2,1
+y

0
v
1,0
= 0, v
2,1
(0) =0,
I. Hashim, M.S.H. Chowdhury / Physics Letters A 372 (2008) 470481 477
(95) v

1,2
T v
1,1
+v
2,1
+3v
2
1,0
v
1,1
=0, v
1,2
(0) =0,
(96) v

2,2
v
1,1
=0, v
2,2
(0) =0,
(97) v

1,3
T v
1,2
+v
2,2
+3v
1,0
v
2
1,1
+3v
2
1,0
v
1,2
=0, v
1,3
(0) =0,
(98) v

2,3
v
1,2
=0, v
2,3
(0) =0.
By the same manipulations as in the previous example, we get
(99) v
1,1
(t ) =T
t
_
0
v
1,0
dt
t
_
0
v
2,0
dt
t
_
0
v
3
1,0
dt =
_
T c
1
c
2
c
3
1
_
t,
(100) v
2,1
(t ) =
t
_
0
v
1,0
dt =c
1
t,
v
1,2
(t ) =T
t
_
0
v
1,1
dt
t
_
0
v
2,1
dt 3
t
_
0
v
2
1,0
v
1,1
dt
(101) =
_
1
2

2
T
2
c
1

1
2

2
T c
2
2
2
T c
3
1

1
2
c
1
+
3
2

2
c
2
1
c
2
+
3
2

2
c
5
1
_
t
2
,
(102) v
2,2
(t ) =
t
_
0
v
1,1
dt =
_
1
2
T c
1

1
2
c
2

1
2
c
3
1
_
t
2
,
v
1,3
(t ) =T
t
_
0
v
1,2
dt
t
_
0
v
2,2
dt 3
t
_
0
_
v
1,0
v
2
1,1
+v
2
1,0
v
1,2
_
dt
=
_
1
6

3
T
3
c
1

5
2

3
c
7
1

1
6

3
T
2
c
2
+3
3
T c
2
1
c
2
+
1
6

2
c
2
+
2
3

2
c
3
1
(103)
1
3

2
T c
1

13
6

3
T
2
c
3
1

3
c
1
c
2
2

7
2

3
c
2
c
4
1
+
9
2

3
T c
5
1
_
t
3
,
(104) v
2,3
(t ) =
t
_
0
v
1,2
dt =
_
1
6

2
T
2
c
1

1
6
T
2
c
2

2
3
T
2
c
3
1

1
6
c
1
+
1
2
c
2
1

2
c
2
+
1
2
c
5
1

2
_
t
3
.
Hence the 4-term approximate series solutions of HPM are
x(t ) =c
1
+
_
T c
1
c
2
c
3
1
_
t +
_
1
2

2
T
2
c
1

1
2
T
2
c
2
2T
2
c
3
1

1
2
c
1
+
3
2
c
2
1

2
c
2
+
3
2
c
5
1

2
_
t
2
+
_
2
3

2
c
3
1
+
1
6

2
c
2

1
3

2
T c
1

5
2

3
c
7
1

7
2

3
c
2
c
4
1
+
9
2

3
T c
5
1

13
6

3
T
2
c
3
1
c
1

3
c
2
2
(105) +3
3
T c
2
1
c
2
+
1
6

3
T
3
c
1

1
6
T
2

3
c
2
_
t
3
,
(106)
y(t ) =c
2
+c
1
t +
_
1
2
T c
1

1
2
c
2

1
2
c
3
1
_
t
2
+
_
1
6

2
T
2
c
1

1
6
T
2
c
2

2
3
T
2
c
3
1

1
6
c
1
+
1
2
c
2
1

2
c
2
+
1
2
c
5
1

2
_
t
3
.
4.2.2. MHPM solution
According to MHPM, we choose the initial approximations as
(107) v
1,0
(t ) =x
0
(t ) =x(t

) =c

1
,
(108) v
2,0
(t ) =y
0
(t ) =y(t

) =c

2
.
Carrying out the steps involved in MHPM gives,
(109) v
1,1
(t ) =T
t
_
t

v
1,0
dt
t
_
t

v
2,0
dt
t
_
t

v
3
1,0
dt =
_
T c

1
c

2
c

1
3
_
(t t

),
478 I. Hashim, M.S.H. Chowdhury / Physics Letters A 372 (2008) 470481
(a) (b)
Fig. 1. Solution of (a) x and (b) y using 4-term HPM, 4-term MHPM (t = 0.001) and RK4 (t = 0.001) related to example 2.
(110) v
2,1
(t ) =
t
_
t

v
1,0
dt =c

1
(t t

),
v
1,2
(t ) =T
t
_
t

v
1,1
dt
t
_
t

v
2,1
dt 3
t
_
t

v
2
1,0
v
1,1
dt
(111) =
_
1
2

2
T
2
c

1

1
2

2
T c

2
2
2
T c

1
3

1
2
c

1
+
3
2

2
c

1
2
c

2
+
3
2

2
c

1
5
_
(t t

)
2
,
(112) v
2,2
(t ) =
t
_
t

v
1,1
dt =
_
1
2
T c

1

1
2
c

2

1
2
c

1
3
_
(t t

)
2
,
v
1,3
(t ) =T
t
_
t

v
1,2
dt
t
_
t

v
2,2
dt 3
t
_
t

_
v
1,0
v
2
1,1
+v
2
1,0
v
1,2
_
dt
=
_
1
6

3
T
3
c

1

5
2

3
c

1
7

1
6

3
T
2
c

2
+3
3
T c

1
2
c

2
+
1
6

2
c

2
+
2
3

2
c

1
3
(113)
1
3

2
T c

1

13
6

3
T
2
c

1
3

3
c

1
c

2
2

7
2

3
c

2
c

1
4
+
9
2

3
T c

1
5
_
(t t

)
3
,
(114) v
2,3
(t ) =
t
_
t

v
1,2
dt =
_
1
6

2
T
2
c

1

1
6
T
2
c

2

2
3
T
2
c

1
3

1
6
c

1
+
1
2
c

1
2

2
c

2
+
1
2
c

1
5

2
_
(t t

)
3
.
To carry out the iterations in every subinterval of equal length, we take the values of the following,
(115) c

1
=x(t

)
4
(t

), c

2
=y(t

)
4
(t

),
where
4
(t ) =

3
i=0
v
1,i
and
4
(t ) =

3
i=0
v
2,i
.
Figs. 1(a) and 1(b) show the solutions for x and y respectively by the 4-term HPM and 4-term MHPM on the time step t =
0.001 and RK4 on the time step t = 0.001 for the case = 40, T = 0.1575, c
1
= 0.45 and c
2
= 0.02025, cf. [24,25]. It is
observed that the 4-term MHPM solutions agree very well with the RK4 solutions for t upto t =5, while the 4-term HPM solutions
are only valid for t 1.
4.3. Example 3
Finally, we shall study the following nonlinear system of ODEs
(116) x

=x(a by),
(117) y

=y(c dx),
I. Hashim, M.S.H. Chowdhury / Physics Letters A 372 (2008) 470481 479
subject to the initial conditions
(118) x(0) =k
1
, y(0) =k
2
,
where x and y are functions of t and a, b, c, d, k
1
, k
2
are known constants.
4.3.1. HPM solution
According to the HPM, we can construct a homotopy of system (116)(117) as follows:
(119) v

1
x

0
+p(x

0
av
1
+bv
1
v
2
) =0,
(120) v

2
y

0
+p(y

0
+cv
2
dv
1
v
2
) =0.
The initial approximations as
(121) v
1,0
(t ) =x
0
(t ) =x(0) =k
1
,
(122) v
2,0
(t ) =y
0
(t ) =y(0) =k
2
.
By the same manipulations as in the previous examples, we have
(123) v
1,1
(t ) =k
1
(a k
2
b)t,
(124) v
2,1
(t ) =k
2
(dk
1
c)t,
(125) v
1,2
(t ) =
1
2
k
1
_
a
2
+2abk
2
k
2
2
b
2
bck
2
+bdk
1
k
2
_
t
2
,
(126) v
2,2
(t ) =
1
2
k
2
_
c
2
2cdk
1
+adk
1
bdk
1
k
2
+d
2
k
1
2
_
t
2
,
v
1,3
(t ) =
1
6
k
1
_
3abck
2
+4abdk
1
k
2
+3b
2
ck
2
2
4b
2
dk
1
k
2
2
3ab
2
k
2
2
(127) +3a
2
bk
2
+bc
2
k
2
2bcdk
1
k
2
+bd
2
k
1
2
k
2
+k
2
3
b
3
a
3
_
t
3
,
v
2,3
(t ) =
1
6
k
2
_
3dc
2
k
1
3cd
2
k
1
2
3acdk
1
+3ad
2
k
1
2
c
3
+4bcdk
1
k
2
(128) 4bd
2
k
1
2
k
2
+a
2
dk
1
2abdk
1
k
2
+b
2
dk
1
k
2
2
+d
3
k
1
3
_
t
3
.
The 4-term approximate solutions for the case a =1, b =1, c =0.1, d =1, k
1
=14 and k
2
=18 are
(129) x(t ) =14 238t +
1358
5
t
2
+
3028697
150
t
3
,
(130) y(t ) =18 +
1251
5
t
40311
100
t
2

20087343
1000
t
3
.
4.3.2. MHPM solution
According to MHPM, we choose the initial approximations as
(131) v
1,0
(t ) =x
0
(t ) =x(t

) =k

1
,
(132) v
2,0
(t ) =y
0
(t ) =y(t

) =k

2
.
Again carrying out the steps involved in MHPM gives
(133) v
1,1
(t ) =k

1
(a k

2
b)(t t

),
(134) v
2,1
(t ) =k

2
(dk

1
c)(t t

),
(135) v
1,2
(t ) =
1
2
k

1
_
a
2
+2abk

2
k

2
2
b
2
bck

2
+bdk

1
k

2
_
(t t

)
2
,
(136) v
2,2
(t ) =
1
2
k

2
_
c
2
2cdk

1
+adk

1
bdk

1
k

2
+d
2
k

1
2
_
(t t

)
2
,
v
1,3
(t ) =
k

1
6
_
3abk

2
c +4abdk

1
k

2
+3b
2
ck

2
2
4b
2
dk

1
k

2
2
3ab
2
k

2
2
a
3
(137) +3a
2
bk

2
+bc
2
k

2
2bcdk

1
k

2
+bd
2
k

1
2
k

2
+k

2
3
b
3
_
(t t

)
3
,
480 I. Hashim, M.S.H. Chowdhury / Physics Letters A 372 (2008) 470481
(a) (b)
Fig. 2. Solutions of (a) x and (b) y using 4-term HPM, 4-term MHPM (t = 0.001) and RK4 (t = 0.001) related to example 3.
v
2,3
(t ) =
k

2
6
_
3c
2
dk

1
3cd
2
k

1
2
3acdk

1
+3ad
2
k

1
2
c
3
+4bcdk

1
k

2
+d
3
k

1
3
(138) +a
2
dk

1
4bd
2
k

1
2
k

2
2abdk

1
k

2
+b
2
dk

1
k

2
2
_
(t t

)
3
.
Figs. 2(a) and 2(b) show the solutions for x and y respectively by the 4-term HPM and 4-term MHPM on the time step t = 0.001
and RK4 on the time step t = 0.001 for the case a = 1, b = 1, c = 0.1, d = 1, k
1
= 14 and k
2
= 18. The same conclusions as in
the second example also apply here.
5. Conclusion
In this Letter, it is shown that the solutions obtained by the standard Homotopy-Perturbation Method (HPM) were not valid for
large time span unless more terms are calculated. We proposed a technique which treated the HPM as an algorithm in a sequence
of intervals (i.e., time step) for getting accurate approximate solutions of linear and nonlinear systems of ODEs. The numerical
comparisons presented demonstrate that the Multistage Homotopy-Perturbation Method (MHPM) is a promising tool for solving
both linear and nonlinear problems.
Acknowledgements
The authors would like to acknowledge the nancial supports received from the Academy of Sciences Malaysia under the
SAGA grant No. P24c (STGL-011-2006), the Malaysian Technical Cooperation Program and the International Islamic University
Chittagong, Bangladesh.
References
[1] S. Olek, SIAM Rev. 36 (1994) 480.
[2] P. Vadasz, S. Olek, Int. J. Heat Mass Transfer 43 (2000) 1715.
[3] I. Hashim, M.S.M. Noorani, R. Ahmad, S.A. Bakar, E.S. Ismail, A.M. Zakaria, Chaos Solitons Fractals 28 (2006) 1149.
[4] M.S.M. Noorani, I. Hashim, R. Ahmad, S.A. Bakar, E.S. Ismail, A.M. Zakaria, Chaos Solitons Fractals 32 (2006) 1296.
[5] O. Abdulaziz, N.F.M. Noor, I. Hashim, M.S.M. Noorani, Chaos Solitons Fractals, doi:10.1016/j.chaos.2006.09.007.
[6] S. Ghosh, A. Roy, D. Roy, Comput. Methods Appl. Mech. Engrg. 196 (2007) 1133.
[7] J. Ruan, Z. Lu, Math. Comput. Modelling, doi:10.1016/j.mcm.2006.12.038.
[8] J.H. He, Non-Perturbative Methods for Strongly Nonlinear Problems, Die Deutsche Bibliothek, Germany, 2006.
[9] J.H. He, Appl. Math. Comput. 135 (2003) 73.
[10] J.H. He, Comput. Methods Appl. Mech. Engrg. 167 (1998) 57.
[11] J.H. He, Int. J. Nonlinear Mech. 35 (2000) 37.
[12] J.H. He, Chaos Solitons Fractals 26 (2005) 695.
[13] J.H. He, Int. J. Nonlinear Sci. Numer. Simul. 6 (2005) 207.
[14] J.H. He, Phys. Lett. A 350 (2006) 87.
[15] J.H. He, Int. J. Mod. Phys. B 20 (2006) 1141.
[16] M.A. Noor, S.T. Mohyud-Din, Math. Comput. Modelling 45 (2007) 954.
[17] A. Rajabi, D.D. Ganji, H. Taherian, Phys. Lett. A 360 (2007) 570.
[18] M. El-Shahed, Int. J. Nonlinear Sci. Numer. Simul. 6 (2005) 63.
[19] D.D. Ganji, A. Sadighi, Int. J. Nonlinear Sci. Numer. Simul. 7 (2006) 411.
[20] D.D. Ganji, A. Sadighi, J. Comput. Appl. Math. 207 (2007) 24.
I. Hashim, M.S.H. Chowdhury / Physics Letters A 372 (2008) 470481 481
[21] M.S.H. Chowdhury, I. Hashim, Phys. Lett. A 365 (2007) 439.
[22] M.S.H. Chowdhury, I. Hashim, O. Abdulaziz, Phys. Lett. A 368 (2007) 251.
[23] M.S.H. Chowdhury, I. Hashim, Phys. Lett. A 368 (2007) 305.
[24] N. Shawagfeh, D. Kaya, Appl. Math. Lett. 17 (2004) 323.
[25] E. Beltrami, Mathematics for Dynamic Modeling, Academic Press, Orlando, 1987.