EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 11, No. 2, 2018, 537-552
ISSN 1307-5543 – www.ejpam.com
Published by New York Business Global

Modifications of the Multistep Optimal Homotopy
Asymptotic Method to some Nonlinear KdV-Equations

M. Fizza1,2, H. Ullah1,2,∗, S. Islam1,2, Q. Shah1,2, F. Chohan1,2, M. Mamat1,2
1 Department of Mathematics, Abdul Wali Khan University Mardan, 23200 Mardan,
Pakistan
2 UET Peshawar, Burraimi University Oman, Sultan Zain al Bidin Malaysia

Abstract. In this article, we have introduced the mathematical theory of Multistep Optimal Homotopy
Asymptotic Method (MOHAM). The proposed method is implemented to different models having system
of partial differential equations (PDEs). The results obtained by the proposed method are compared with
Homotopy Analysis Method (HAM) and closed form solutions. The comparisons of these results show
that (MOHAM) is simpler in applicability, effective, explicit, control the convergence through optimal
constants, involve less computational work. The (MOHAM) is independent of the assumption of initial
conditions and small parameters like Homotopy Perturbation Method (HPM), Homotopy Analysis Method
(HAM), Variational Iteration Method (VIM), Adomian Decomposition Method (ADM) and Perturbation
Method (PM).

2010 Mathematics Subject Classifications: 35C05, 35Q53

Key Words and Phrases: Multistep Optimal Homotopy Asymptotic Method (MOHAM), closed form
solutions, KdV equations.

1. Introduction

The nonlinear physical phenomena can be understood by the effort of finding the exact so-
lution. The exact solutions of some boundary value problems (BVPs) may be found by using
transformation based techniques like invariance group analysis method [1], Lie infinitesimal crite-
rion [2], the symbolic computation [3] and Backlund transformation [4]. These methods reduced
the complex equations into simple equations by using the transformation. The exact solutions
of all the nonlinear problems are either not available or difficult to find because (PDEs) have
infinitely many solutions. So the attentions of the researchers are attracted to develop the ap-
proximation tools for the nonlinear (BVPs) like Variational Iteration Method (VIM) [5], Adomian
Decomposition Method (ADM) [6], Differential Transform Method (DTM) [7], Homotopy Pertur-
bation Method (HPM) [8] and Perturbation Method (PM) [9–11] have been used for the solution
of nonlinear (PDEs). These methods contain a small parameter which cannot be found easily.
The homotopy was combined with perturbation techniques such as Homotopy Analysis Method
(HAM) [12] and Homotopy Perturbation Method (HPM) [8]. For these methods an initial solu-
tion is needed to assume.
Marinca et al. introduced (OHAM) [14–16] for the solution of nonlinear problems which made
the perturbation methods independent of the assumption of small parameters and huge compu-
tational work.
The motivation of this paper is to formulate and implement the (MOHAM) for the solution of
system of (PDEs). In [19–23], (OHAM) has been proved to be valuable for obtaining approximate
solutions of various boundary value problems. Here we have proved that (MOHAM) is useful

∗Corresponding author.

Email addresses: hakeemullah1@gmail.com (H. Ullah), ahmadjan524@gmail.com (M. Fiza), saeed.sn@gmail.com
(S. Islam), qayyumshah@uetpeshawar.edu.pk (Q. Shah),
farkhandayusaf@yahoo.com (F.Chohan), musmat567@gmail.com (M.Mamat)

http://www.ejpam.com 537 c© 2018 EJPAM All rights reserved.



H. Ullah et al. / Eur. J. Pure Appl. Math, 11 (2) (2018), 537-552 538

and reliable for the complex nonlinear (PDEs), showing its validity and great potential for the
solution of transient physical phenomenon in science and engineering.
In the succeeding section, the basic idea of (MOHAM) is formulated. The effectiveness and
efficiency of (MOHAM) is shown in Section 3.

2. Fundamental Mathematical Theory of (MOHAM) for a system of coupled
PDEs

Consider a system of n partial differential equations of the form

A1

(
f1(ζ, t), f2(ζ, t), ..., fn(ζ, t)

)
+ s1(ζ, t) = 0,

A2

(
f1(ζ, t), f2(ζ, t), ..., fn(ζ, t)

)
+ s2(ζ, t) = 0,

.

.

.

An

(
f1(ζ, t), f2(ζ, t), ..., fn(ζ, t)

)
+ sn(ζ, t) = 0, ζ ∈ Ω

(1)

with boundary conditions 

B1

(
f1(ζ, t),

df1(ζ, t)

dζ

)
= 0,

B2

(
f2(ζ, t),

df2(ζ, t)

dζ

)
= 0,

.

.

.

Bn

(
fn(ζ, t),

dfn(ζ, t)

dζ

)
= 0, ζ ∈ Γ

(2)

where A1, A2, ..., An are differential operators, f1(ζ, t), f2(ζ, t), ..., fn(ζ, t) are unknown func-
tions, ζ and t denote spatial variables, respectively, Γ is the boundaries of Ωand s1(ζ, t), s2(ζ, t), ..., sn(ζ, t)
are known analytic functions.
A1, A2, ..., An can be divided into two parts:

A1 = L1 +N1,

A2 = L2 +N2,

.

.

.

An = Ln +Nn

(3)

L1, L2, ..., Ln contain the linear parts while N1, N2, ..., Nn contains the nonlinear parts of the
system of partial differential equations.



H. Ullah et al. / Eur. J. Pure Appl. Math, 11 (2) (2018), 537-552 539

According to (OHAM), one can construct a family of equations

(1− r)
[
L1

(
ϕ(ζ, t, r)

)
+ s1(ζ, t)

]
= H1(r)

[
L1

(
ϕ(ζ, t, r)

)
+N1

(
ϕ(ζ, t, r)

)
+ s1(ζ, t)

]
= 0,

(1− r)
[
L2

(
ϕ(ζ, t, r)

)
+ s2(ζ, t)

]
= H2(r)

[
L2

(
ϕ(ζ, t, r)

)
+N2

(
ϕ(ζ, t, r)

)
+ s2(ζ, t)

]
= 0,

.

.

.

(1− r)
[
Ln

(
ϕ(ζ, t, r)

)
+ sn(ζ, t)

]
= Hn(r)

[
Ln

(
ϕ(ζ, t, r)

)
+Nn

(
ϕ(ζ, t, r)

)
+ sn(ζ, t)

]
= 0

(4)
where the auxiliary functions H1((ζ, t, r)), H2((ζ, t, r)), ...,Hn((ζ, t, r)) are nonzero for r 6= 0.

(4) is called optimal homotopy equations. Clearly, we have

r = 0⇒ L1

(
ϕ(ζ, t, 0)

)
+ s1(ζ, t) = 0,

r = 0⇒ L2

(
ϕ(ζ, t, 0)

)
+ s2(ζ, t) = 0,

.

.

.

r = 0⇒ Ln

(
ϕ(ζ, t, 0)

)
+ sn(ζ, t) = 0

(5)

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

r = 1⇒ L1

(
ϕ(ζ, t, 1)

)
+ s1(ζ, t) = 0,

r = 1⇒ L2

(
ϕ(ζ, t, 1)

)
+ s2(ζ, t) = 0,

.

.

.

r = 1⇒ Ln

(
ϕ(ζ, t, 1)

)
+ sn(ζ, t) = 0

(6)

Obviously, when r = 0 and r = 1, we obtain

ϕ1(ζ, t, 0) = (f1)0(ζ, t),

ϕ2(ζ, t, 0) = (f2)0(ζ, t),

.

.

.

ϕn(ζ, t, 0) = (fn)0(ζ, t)

(7)

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

ϕ1(ζ, t, 1) = f1(ζ, t),

ϕ2(ζ, t, 1) = f2(ζ, t),

.

.

.

ϕn(ζ, t, 1) = fn(ζ, t)

(8)

respectively.
When r varies from 0 to 1, the solution ϕ1(ζ, t, r), ϕ2(ζ, t, r), ..., ϕn(ζ, t, r) approaches from



H. Ullah et al. / Eur. J. Pure Appl. Math, 11 (2) (2018), 537-552 540

(f1)0(ζ, t), (f2)0(ζ, t), ..., (fn)0(ζ, t) to f1(ζ, t), f2(ζ, t), ..., f1(ζ, t)

L1

(
(f1)0)(ζ, t)

)
+ s1(ζ, t) = 0, B1

(
(f1)0(ζ, t),

d(f1)0(ζ, t)

dζ

)
= 0,

L2

(
(f2)0)(ζ, t)

)
+ s2(ζ, t) = 0, B2

(
(f2)0(ζ, t),

d(f2)0(ζ, t)

dζ

)
= 0,

.

.

.

Ln

(
(fn)0)(ζ, t)

)
+ sn(ζ, t) = 0, Bn

(
(fn)0(ζ, t),

d(fn)0(ζ, t)

dζ

)
= 0

(9)

We choose auxiliary functions H1(r), H2(r), ...,Hn(r) in the form

H1(r) = rC11 + r2C12 + ...+ rnC1n,

H2(r) = rC21 + r2C22 + ...+ rnC2n,

.

.

.

Hn(r) = rCn1 + r2Cn2 + ...+ rnCnn

(10)

To get the approximate solutions, we expand ϕ1(ζ, t, r, C1i), ϕ2(ζ, t, r, C2i), ..., ϕn(ζ, t, r, Cni)
by Taylor’s series about γ in the following manner,

ϕ1(ζ, t, r, C1i) = (f1)0(ζ, t) +
∑
k≥1

(f1)k(ζ, t, C1i)r
k,

ϕ2(ζ, t, r, C2i) = (f2)0(ζ, t) +
∑
k≥1

(f2)k(ζ, t, C2i)r
k,

.

.

.

ϕn(ζ, t, r, Cni) = (fn)0(ζ, t) +
∑
k≥1

(fn)k(ζ, t, Cni)r
k

(11)

where k = 1, 2, .... Now substituting Equations (10-11) into Equation (4) and equating the
coefficient of like powers of r, we obtain Zeroth order system, given by Equation (9), the first and
second order systems given by Equations (12-14) respectively and the general governing equations
are given by Equation (15)

L1

(
(f1)1(ζ, t)

)
− L1

(
(f1)0(ζ, t)

)
= C11

(
L1

(
(f1)0(ζ, t)

)
+N1

(
(f1)0(ζ, t)

))
, B1

(
(f1)1(ζ, t),

d(f1)1(ζ, t)

dζ

)
,

L2

(
(f2)1(ζ, t)

)
− L2

(
(f2)0(ζ, t)

)
= C21

(
L2

(
(f2)0(ζ, t)

)
+N2

(
(f2)0(ζ, t)

))
, B2

(
(f2)1(ζ, t),

d(f2)1(ζ, t)

dζ

)
,

.

.

.

Ln

(
(fn)1(ζ, t)

)
− Ln

(
(fn)0(ζ, t)

)
= Cn1

(
Ln

(
(fn)0(ζ, t)

)
+Nn

(
(fn)0(ζ, t)

))
, Bn

(
(fn)1(ζ, t),

d(fn)1(ζ, t)

dζ

)
(12)



H. Ullah et al. / Eur. J. Pure Appl. Math, 11 (2) (2018), 537-552 541

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

L1

(
(f1)2(ζ, t)

)
− L1

(
(f1)1(ζ, t)

)
= C11

(
L1

(
(f1)1(ζ, t)

))
+N1

(
(f1)1(ζ, t), (f1)0(ζ, t)

)
+ C12

(
L1

(
(f1)0(ζ, t)

)
+N1

(
(f1)0(ζ, t)

))
, B1

(
(f1)2(ζ, t),

d(f1)2(ζ, t)

dζ

)
= 0,

L2

(
(f2)2(ζ, t)

)
− L2

(
(f2)1(ζ, t)

)
= C21

(
L2

(
(f2)1(ζ, t)

))
+N2

(
(f2)1(ζ, t), (f2)0(ζ, t)

)
+ C22

(
L2

(
(f2)0(ζ, t)

)
+N2

(
(f2)0(ζ, t)

))
, B2

(
(f2)2(ζ, t),

d(f2)2(ζ, t)

dζ

)
= 0,

.

.

.

Ln

(
(fn)n(ζ, t)

)
− Ln

(
(fn)1(ζ, t)

)
= Cn1

(
Ln

(
(fn)1(ζ, t)

))
+Nn

(
(fn)1(ζ, t), (fn)0(ζ, t)

)
+ Cn2

(
Ln

(
(fn)0(ζ, t)

)
+Nn

(
(fn)0(ζ, t)

))
, Bn

(
(fn)2(ζ, t),

d(fn)2(ζ, t)

dζ

)
= 0

(13)

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

L1

(
(f1)3(ζ, t)

)
− L1

(
(f1)2(ζ, t)

)
= C11

(
L1

(
(f1)2(ζ, t)

)
+N1

(
(f1)2(ζ, t), (f1)1(ζ, t), (f1)0(ζ, t)

))
+ C12

(
L1

(
(f1)1(ζ, t)

)
+N1

(
(f1)1(ζ, t), (f1)0(ζ, t)

))
+ C13

(
L1

(
(f1)0(ζ, t)

)
+N1

(
(f1)0(ζ, t)

))
, B1

(
(f1)3(ζ, t),

d(f1)3(ζ, t)

dζ

)
= 0,

L2

(
(f2)3(ζ, t)

)
− L2

(
(f2)2(ζ, t)

)
= C21

(
L2

(
(f2)2(ζ, t)

)
+N2

(
(f2)2(ζ, t), (f2)1(ζ, t), (f2)0(ζ, t)

))
+ C22

(
L2

(
(f2)1(ζ, t)

)
+N2

(
(f2)1(ζ, t), (f2)0(ζ, t)

))
+ C23

(
L2

(
(f2)0(ζ, t)

)
+N2

(
(f2)0(ζ, t)

))
, B2

(
(f2)3(ζ, t),

d(f2)3(ζ, t)

dζ

)
= 0,

.

.

.

Ln

(
(fn)3(ζ, t)

)
− Ln

(
(fn)2(ζ, t)

)
= Cn1

(
Ln

(
(fn)2(ζ, t)

)
+Nn

(
(fn)2(ζ, t), (fn)1(ζ, t), (fn)0(ζ, t)

))
+ Cn2

(
Ln

(
(fn)1(ζ, t)

)
+Nn

(
(fn)1(ζ, t), (fn)0(ζ, t)

))
+ Cn3

(
Ln

(
(fn)0(ζ, t)

)
+Nn

(
(fn)0(ζ, t)

))
, Bn

(
(fn)3(ζ, t),

d(fn)3(ζ, t)

dζ

)
= 0

(14)



H. Ullah et al. / Eur. J. Pure Appl. Math, 11 (2) (2018), 537-552 542

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

L1

(
(f1)k(ζ, t)

)
− L1

(
(f1)k−1(ζ, t)

)
=

k∑
i=1

C1i

[
L1

(
(f1)k−1(ζ, t)

)
+N1

(
(f1)k−1(ζ, t), (f1)k−2(ζ, t), ..., (f1)0(ζ, t)

)]
,

k = 2, 3, ..., B1

(
(f1)k(ζ, t),

d(f1)k(ζ, t)

dζ

)
= 0,

L2

(
(f2)k(ζ, t)

)
− L2

(
(f2)k−1(ζ, t)

)
=

k∑
i=1

C2i

[
L2

(
(f2)k−1(ζ, t)

)
+N2

(
(f2)k−1(ζ, t), (f2)k−2(ζ, t), ..., (f2)0(ζ, t)

)]
,

k = 2, 3, ..., B2

(
(f2)k(ζ, t),

d(f2)k(ζ, t)

dζ

)
= 0,

.

.

.

Ln

(
(fn)k(ζ, t)

)
− Ln

(
(fn)k−1(ζ, t)

)
=

k∑
i=1

Cni

[
Ln

(
(fn)k−1(ζ, t)

)
+Nn

(
(fn)k−1(ζ, t), (fn)k−2(ζ, t), ..., (fn)0(ζ, t)

)]
,

k = 2, 3, ..., Bn

(
(fn)k(ζ, t),

d(fn)k(ζ, t)

dζ

)
= 0

(15)
It has been observed that the convergence of the series Equations (11) depends upon the auxiliary
constants C11, C12, ..., C21, C22, ..., Cn1, Cn2, ....
If it is convergent at , one has

ϕ1(ζ, t, C1i) = (f1)0)(ζ, t) +
∑
k≥1

(f1)k(ζ, t, C1i), i = 1, 2, ...,m,

ϕ2(ζ, t, C2i) = (f2)0)(ζ, t) +
∑
k≥1

(f2)k(ζ, t, C1i), i = 1, 2, ...,m,

.

.

.

ϕn(ζ, t, Cni) = (fn)0)(ζ, t) +
∑
k≥1

(fn)k(ζ, t, Cni), i = 1, 2, ...,m

(16)

Generally speaking the solution of Equations (1) can be approximately written as

fm1 (ζ, t, C1i) = (f1)0)(ζ, t) +
m∑
k≥1

(f1)k(ζ, t, C1i), i = 1, 2, ...,

f2M(ζ, t, C2i) = (f2)0)(ζ, t) +

m∑
k≥1

(f2)k(ζ, t, C1i), i = 1, 2, ...,

.

.

.

fmn (ζ, t, Cni) = (fn)0)(ζ, t) +

m∑
k≥1

(fn)k(ζ, t, Cni), i = 1, 2, ...

(17)

where m is the order of approximation. Substituting Equations (17) into Equations (1), it results



H. Ullah et al. / Eur. J. Pure Appl. Math, 11 (2) (2018), 537-552 543

the following expression for residuals

R1(ζ, t, C1i) = L1

(
fm1 (ζ, t, C1i)

)
+ s1(ζ, t) +N1

(
fm1 (ζ, t, C1i)

)
,

R2(ζ, t, C2i) = L2

(
fm2 (ζ, t, C2i)

)
+ s2(ζ, t) +N2

(
fm2 (ζ, t, C2i)

)
.

.

.

Rn(ζ, t, Cni) = Ln

(
fmn (ζ, t, Cni)

)
+ sn(ζ, t) +Nn

(
fmn (ζ, t, Cni)

)
(18)

IfR1(ζ, t, C1i) = 0, R2(ζ, t, C2i) = 0, ..., Rn(ζ, t, Cni) = 0, then ϕ1(ζ, t, C1i), ϕ2(ζ, t, C2i), ..., ϕn(ζ, t, Cni)
will be the exact solutions of the problem. Generally it doesnt happen, especially in nonlinear
problems.
For the computation of auxiliary constants, C1i, C2i, ..., Cni, i = 1, 2, ...,m, there are different
methods like Galerkins Method, Ritz Method, Least Squares Method and Collocation Method.
One can apply the Method of Least Squares as under

J1(C1i) =

∫
Γ

∫
Ω
R2

1(ζ, t, C1i) dζ dt,

J2(C2i) =

∫
Γ

∫
Ω
R2

2(ζ, t, C2i) dζ dt,

.

.

.

Jn(Cni) =

∫
Γ

∫
Ω
R2
n(ζ, t, Cni) dζ dt

(19)

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

∂J1

∂C11
=

∂J1

∂C12
= ... =

∂J1

∂C1m
,

∂J2

∂C21
=

∂J2

∂C22
= ... =

∂J2

∂C2m
,

.

.

.

∂Jn
∂Cn1

=
∂Jn
∂Cn2

= ... =
∂Jn
∂Cnm

(20)

The m th order approximate solution can be obtained by these constants. The more general
auxiliary functions H1(r), H2(r), ...,Hn(r) are useful for convergence, which depends upon con-
stants C11, C12, ..., Cn1, Cn2, ..., can be optimally identified by Equations (19) and is useful in
error minimization.
The same procedure is repeating for the next iteration and so on.
The (MOHAM) solution is given by

f̃1(ζ) =


f1(ζ) , z0 ≤ t ≤ z1

.........
fN (ζ) , zN−1 ≤ t ≤ T

f̃2(ζ) =


f2(ζ) , z0 ≤ t ≤ z1

.........
fN (ζ) , zN−1 ≤ t ≤ T



H. Ullah et al. / Eur. J. Pure Appl. Math, 11 (2) (2018), 537-552 544

.

.

.

f̃n(ζ) =


fn(ζ) , z0 ≤ t ≤ z1

.........
fN (ζ) , zN−1 ≤ t ≤ T

3. Implementation of the MOHAM formulation of a system of PDEs

Model(1): Application of (MOHAM) to a system of three (PDEs) KdV equations
of the form 

∂u(ζ, t)

∂t
= 0.5

∂3u(ζ, t)

∂ζ3
− 3u(ζ, t)

∂u(ζ, t)

∂ζ
+ 3

∂

∂ζ

(
v(ζ, t)w(ζ, t)

)
,

∂v(ζ, t)

∂t
= −∂

3v(ζ, t)

∂3ζ3
+ 3u(ζ, t)

∂v(ζ, t)

∂ζ
,

∂w(ζ, t)

∂t
= −∂

3w(ζ, t)

∂ζ3
+ 3u(ζ, t)

∂w(ζ, t)

∂t

(21)

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

u(ζ, 0) =
1

3
(β − 8k2) + 4k2 tanh2(kζ),

v(ζ, 0) =
−4k2

(
3k2c0 − 2βc2 + 4k2c2

)
3c2

2

+
4k2

c2
tanh2(kζ),

w(ζ, 0) = c0 + c2 tanh2(kζ)

(22)

The closed form solution of Equations (21) is given by [24]

u(ζ, t) =
1

3
(β − 8k2) + 4k2 tanh2(k(ζ + βt)),

v(ζ, t) =
−4k2

(
3k2c0 − 2βc2 + 4k2c2

)
3c2

2

+
4k2

c2
tanh2(k(ζ + βt)),

w(ζ, 0) = c0 + c2 tanh2(k(ζ + βt))

(23)

Applying the formulation of Extended (OHAM) technique discussed in section 2.
We consider 

u = u0 + γu1 + γ2u2, v = v0 + γv1 + γ2v2,

w = w0 + γw1 + γ2w2, H1(γ) = γC11 + γ2C12,

H2(γ) = γc21 + γ2C22, H3(γ) = γC31 + γ2C32

(24)

Zeroth Order System:

2
∂u0

∂t
= 0,

∂v0

∂t
= 0,

∂w0

∂t
= 0, (25)

with initial conditions

u0(ζ, 0) =
1

3
(β − 8k2) + 4k2 tanh2(kζ),

v0(ζ, 0) =
−4k2

(
3k2c0 − 2βc2 + 4k2c2

)
3c2

2

+
4k2

c2
tanh2(kζ),

w0(ζ, 0) = c0 + c2 tanh2(kζ)

(26)



H. Ullah et al. / Eur. J. Pure Appl. Math, 11 (2) (2018), 537-552 545

Its solution is 

u0(ζ, t) =
1

3

(
β − 8k2 + 12k2 tanh2(kζ)

)
,

v0(ζ, t) =
−4
(

3k2c0 + 4k2c2 − 2k2βc2 − 3k2c2 tanh2(kζ)
)

3c2
2

,

w0(ζ, t) = c0 + c2 tanh2(kζ)

(27)

First Order System:

2
∂u1(ζ, t)

∂t
= 2(1 + C11)

∂u0

∂t
+ 6C11

(
u0
∂u0

∂ζ
− w0

∂v0

∂ζ
− v0

∂w0

∂ζ

)
− C11

∂3u0

∂ζ3
,

∂v1(ζ, t)

∂t
= (1 + C21)

∂v0

∂t
− 3C21u0

∂w0

∂ζ
− C21

∂3w0

∂ζ3
,

∂w1(ζ, t)

∂t
= (1 + C31)

∂w0

∂t
− 3C31w0

∂w0

∂ζ
+ C31

∂3w0

∂ζ3

(28)

with

u1(ζ, 0) = 0, v1(ζ, 0) = 0, w1(ζ, 0) = 0 (29)

Its solution is

u1(ζ, t, C11) =
8tC11

c2

[
− 3k2sech2(kζ)c0 tanh(kζ) + 3k5sech2(kζ)c0 tanh(kζ)− 4k2sech2(kζ)c2 tanh(kζ)

− k3βsech2(kζ)c2 tanh(kζ) + 4k5sech4(kζ)c2 tanh(kζ)− 6k3sech2(kζ)c2 tanh2(kζ) + 10k3sech2(kζ)c2 tanh3(kζ)

]
,

v1(ζ, t, C21) = −8tC21

c2

[
− 8k5sech2(kζ) tanh(kζ) + k3βsech2(kζ) tanh(kζ) + 8k5sech4(kζ) tanh(kζ)

+ 8k5sech5(kζ) tanh3(kζ)

]
,

w1(ζ, t, C31) = −2tC31

[
− 8k3sech2(kζ)c2 tanh(kζ) + kβsech2(kζ)c2 tanh(kζ) + 8k3sech4(kζ)c2 tanh(kζ)

+ 8k3sech2(kζ)c2 tanh3(kζ)

]
(30)

The approximate solution is obtained as
u(ζ, t, C11) = u0(ζ, t) + u1(ζ, t, C11),

v(ζ, t, C21) = v0(ζ, t) + v1(ζ, t, C21),

w(ζ, t, C31) = w0(ζ, t) + w1(ζ, t, C31)

(31)

For the computation of the constants C11, C21, and C31 using (31) in (21) and applying the
technique mentioned in (18 - 20) by taking c0 = 1.5, c2 = 0.1, k = 0.1, β = 1.5 and also
repeating the same procedure for next iteration, We get the approximate solution as

u(ζ, t) =

{
1
3 [1.42 + 0.12 tanh2(0.1ζ)] + 80t[4.210289× 10−14sech2(0.1ζ) tanh(0.1ζ)] , 0 ≤ t ≤ 0.5
80t[−3.65449× 10−17sech4(0.1ζ) tanh(0.1ζ) + 5.39038× 10−15sech2(0.1ζ) tanh3(0.1ζ)] , 0.5 ≤ t ≤ 1



H. Ullah et al. / Eur. J. Pure Appl. Math, 11 (2) (2018), 537-552 546

v(ζ, t) =

{
−133.333(−0.00051− 0.001 tanh2(0.1ζ)) , 0 ≤ t ≤ 0.5
−133.333(−0.00051− 0.002 tanh2(0.1ζ)) , 0.5 ≤ t ≤ 1

w(ζ, t) =

{
1.5 + 0.1 tanh2(0.1ζ)− 2t[−2.65863× 10−21sech2(0.1ζ) tanh(0.1ζ)] , 0 ≤ t ≤ 0.5
−2t[−1.49782× 10−22sech4(0.1ζ) tanh(0.1ζ)− 1.49782× 10−22sech2(0.1ζ) tanh2(0.1ζ)] , 0.5 ≤ t ≤ 1

for C11 = −9.136232469244701 × 10−12, C12 = C21 = C22 = 0, C31 = −1.8722771193005665 ×
10−19, C32 = 0

Table 1: Comparison of (MOHAM), HAM and Closed form solutions for u(ζ, t) at t = 1

ζ (MOHAM) Solution Closed Form Solution (HAM) Solution

0 0.473333 0.47422 0.473333

20 0.510507 0.51122 0.51030

40 0.51328 0.513294 0.51324

60 0.513332 0.513333 0.513310

80 0.513333 0.513333 0.513333

100 0.513333 0.513333 0.513333

Table 2: Comparison of (MOHAM), HAM and Closed form solutions for v(ζ, t) at t = 1

ζ (MOHAM) Solution Closed Form Solution (HAM) Solution

0 0.334667 0.343533 0.334661

20 0.706406 0.713534 0.706402

40 0.73413 0.734269 0.73411

60 0.734657 0.734659 0.734654

80 0.734666 0.734667 0.734665

100 0.734667 0.734667 0.734667

Table 3: Comparison of (MOHAM), HAM and Closed form solutions for w(ζ, t) at t = 1

ζ (MOHAM) Solution Closed Form Solution (HAM) Solution

0 1.5 1.50222 1.5

20 1.59393 1.59472 1.59293

40 1.59987 1.5999 1.59985

60 1.6 1.6 1.6

80 1.6 1.6 1.6

100 1.6 1.6 1.6



H. Ullah et al. / Eur. J. Pure Appl. Math, 11 (2) (2018), 537-552 547

Table 4: Absolute error of (MOHAM) solution u(ζ, t) corresponding to the Closed form Solution

ζ t = 1 t = 0.5 t = 0.1 t = 0.01

0 3.27717×10−2 1.61366×10−2 8.8667×10−4 8.99865 ×10−6

20 2.6804×10−3 2.17746×10−3 7.128×10−4 8.06039×10−5

40 5.09658×10−5 4.16635×10−5 1.38951×10−5 1.58421×10−6

60 9.34118×10−7 7.63709×10−7 2.54789×10−7 2.90535×10−8

80 1.71092×10−8 1.3988×10−8 4.66673×10−9 5.32147×10−10

100 3.1336×10−10 2.562×10−10 8.54742×10−11 9.74665×10−12

Table 5: Absolute error of (MOHAM) solution v(ζ, t) corresponding to the Closed form Solution

ζ t = 1 t = 0.5 t = 0.1 t = 0.01

0 0.327717 0.161366 8.8667×10−3 8.99865×10−5

20 2.6804×10−2 2.17746×10−2 7.128×10−3 8.06039×10−4

40 5.09658×10−4 4.16635×10−4 1.38951×10−4 1.58421×10−5

60 9.34118×10−6 7.63709×10−6 2.54789×10−6 2.90535×10−7

80 1.71092×10−7 1.3988×10−7 4.66673×10−8 5.32146×10−9

100 3.13366×10−9 2.562×10−9 8.54741×10−10 9.7466×10−11

Model(2): Application of the (MOHAM) to Kdv Equation of the form:
∂u(ζ, t)

∂t
= 0.5uζζζ − 3u2uζ + 1.5vζζ + 3uvζ + 3uζv − 3λuζ ,

∂v(ζ, t)

∂t
= −vζζζ − 3vvζ − 3uζvζ + 3u2vζ + 3λvζ

(32)

with initial conditions {
u(ζ, 0) = γ tanh(γζ),

v(ζ, 0) = 0.5(4γ2 + λ)− 2γ2 tanh2(γζ)
(33)

The Closed form solution of the problem is given as [25]
u(ζ, t) = γ tanh

(
γζt(γ2 + 1.5λ)

)
,

v(ζ, t) = 0.5(4γ2 + λ)− 2γ2 tanh2

(
γ
(
ζ − (γ2 + 1.5λ)t

)) (34)

Zeroth Order System:

2
∂u0

∂t
= 0,

∂v0

∂t
= 0, (35)

with initial conditions {
u0(ζ, 0) = γ tanh(γζ,

v0(ζ, 0) = 0.5(4γ2 + λ)− 2γ2 tanh2(γζ)
(36)

Its solution is {
u0(ζ, t) = γ tanh(γζ,

v0(ζ, t) = 0.5(4γ2 + λ)− 2γ2 tanh2(γζ)
(37)

First Order System:
2
∂u1(ζ, t)

∂t
= 2(1 + C11)

∂u0

∂t
+ 6C11(λ− u2

0)
∂u0

∂ζ
+ 6C11u0

∂v0

∂ζ
+ 3C11

∂2v0

∂ζ2
+ C11

∂3u0

∂ζ3
,

∂v1

∂t
= (1 + C21)

∂v0

∂t
− 3C11(λ+ u2

0)
∂v0

∂ζ
− 3(v0 +

∂u0

∂ζ
)
∂v0

∂ζ
− C21

∂3v0

∂ζ3

(38)



H. Ullah et al. / Eur. J. Pure Appl. Math, 11 (2) (2018), 537-552 548

Table 6: Absolute error of (MOHAM) solution w(ζ, t) corresponding to the Closed form Solution

ζ t = 1 t = 0.5 t = 0.1 t = 0.01

0 8.19293×10−2 4.03414×10−2 2.21668×10−3 2.24966×10−5

20 6.70099×10−3 5.44365×10−3 1.782×10−3 2.0151×10−4

40 1.27415×10−4 1.04159×10−4 3.47377×10−5 3.96053×10−6

60 2.33529×10−6 1.90927×10−6 6.36974×10−7 7.26338×10−8

80 4.27729×10−8 3.49701×10−8 1.16668×10−8 1.33037×10−9

100 7.83414×10−10 6.40499×10−10 2.13686×10−10 2.43665×10−11

(a) 2D - Residual of u(ζ, t) at t = 1

0 20 40 60 80

0.00

0.05

0.10

0.15

0.20

0.25

0.30

x

R
e
s
id

u
a
l

u

Residual u

(b) 2D - Residual of v(ζ, t) at t = 1

0 20 40 60 80

0.00

0.01

0.02

0.03

0.04

x

R
e
s
id

u
a
l

v

Residual v

(c) 2D - Residual of w(ζ, t) at t = 1

0 20 40 60 80

0.000

0.002

0.004

0.006

0.008

0.010

x

R
e
s
id

u
a
l

w

Residual w

Figure 1: 2D - Residual

with initial conditions

u1(ζ, 0) = 0, v1(ζ, 0) = 0, (39)

Its solution is 
u1(ζ, t, C11) =

tγsech2(γζ)

2

[
(2γ2 − 3λ)C11)− 6λC11

]
,

v1(ζ, t, C21) = 2tγ3C21

[
− 2γ2 + 3λ

]
sech2(γζ) tanh(γζ)

(40)

The approximate solution is obtained as{
u1(ζ, t, C11) = u0(ζ, t) + u1(ζ, t, C11),

v1(ζ, t, C21) = v0(ζ, t) + v1(ζ, t, C21)
(41)

Repeating the same step for next iteration we obtained the approximate solution as

u(ζ, t) =

{
−0.0002304230sech4(0.1ζ) + 0.1 tanh(0.1ζ) , 0 ≤ t ≤ 0.5
tsech2(0.1ζ)(0.00150104− 0.0000230423 tanh2(0.1ζ)) , 0.5 ≤ t ≤ 1



H. Ullah et al. / Eur. J. Pure Appl. Math, 11 (2) (2018), 537-552 549

v(ζ, t) =

{
0.52 + 0.0000605366tsech4(0.1ζ) tanh(0.1ζ)− 0.02 tanh2(0.1ζ) , 0 ≤ t ≤ 0.5
tsech2(0.1ζ)(0.0174345 + 0.00005366 tanh2(0.1ζ)) , 0.5 ≤ t ≤ 1

for C11 = −0.3291761636279658, C12 = −0.00007616362, C21 = 3.0068307906034297, C22 =
0.020831245, γ = 0.1, λ = 1.

Table 7: Absolute error of (MOHAM) solution u(ζ, t) corresponding to the Exact Solution

ζ t = 1 t = 0.5 t = 0.1 t = 0.01

0 1.64508×10−2 8.144×10−3 4.49177×10− 4.55951×10−9

20 2.17387×10−3 1.90041×10−3 7.54522×10−4 8.11496×10−5

40 2.73879×10−4 2.65752×10−5 1.47276×10−5 1.59572×10−6

60 5.14463×10−6 4.83366×10−7 2.70061×10−7 2.92649×10−8

80 9.42707×10−8 8.85201×10−9 4.94645×10−9 5.36017×10−10

100 1.72664×10−9 1.6213×10−10 9.05975×10−11 9.81748×10−12

Table 8: Absolute error of (MOHAM) solution v(ζ, t) corresponding to the Exact Solution

ζ t = 1 t = 0.5 t = 0.1 t = 0.01

0 1.64508×10−2 8.144×10−3 4.49177×10− 4.55951×10−9

20 2.17387×10−3 1.90041×10−3 7.54522×10−4 8.11496×10−5

40 2.73879×10−4 2.65752×10−5 1.47276×10−5 1.59572×10−6

60 5.14463×10−6 4.83366×10−7 2.70061×10−7 2.92649×10−8

80 9.42707×10−8 8.85201×10−9 4.94645×10−9 5.36017×10−10

100 1.72664×10−9 1.6213×10−10 9.05975×10−11 9.81748×10−12

(a) Comparison of MOHAM and Exact solutions at t =
1

0 20 40 60 80 100

0.086

0.088

0.090

0.092

0.094

0.096

0.098

0.100

x

u

Exact

HPM

MOHAM

(b) Comparisons of MOHAM and Exact solutions at
t = 1

0 20 40 60 80 100

0.500

0.501

0.502

0.503

0.504

0.505

0.506

x

v

Exact sol

HPM

MOHAM

Figure 2:



REFERENCES 550

(a) 2D - Residual of u(ζ, t) at t = 1

0 20 40 60 80 100

-0.001

0.000

0.001

0.002

0.003

0.004

0.005

x

R
e
s
id

u
a
l

u

Residual u

(b) 2D - Residual of u(ζ, t) at t = 1

0 20 40 60 80

-0.0001

0.0000

0.0001

0.0002

x

R
e
s
id

u
a
l

v

Residual v

Figure 3: 2D - Residual

4. Results and Discussions

The mathematical theory of (MOHAM) provides highly accurate solutions for the system of
BVP presented in section 3. We have used Mathematica 7 for our computational work. In Table
1-2, the (MOHAM) results are compared with closed form and (HAM) solutions at t = 1 for the
generalized Hirota-Satsuma coupled KdV equation. The absolute errors at different values of t
are given in Tables 1-6. The residuals are plotted in Figure 1. While for the (MKDV) equation,
the Figure 2 shows the accuracy of (MOHAM) for t = 1. The Tables 7-8 presents the absolute
errors at different values of t. The residuals of MKDV equations are plotted in Figure 3 at t = 1.
From these Tables and Figure. it is evident that the (MOHAM) results are more accurate than
(HAM) results and nearly identical to closed form solutions not only for small values of t but for
large values of t. Here the results are very consistent with the decreasing time and by increasing
displacement.

5. Conclusion

In this paper, we have seen the effectiveness of (MOHAM) different models of system of PDEs.
The (MOHAM) is simpler in applicability, more convenient to control convergence and involved
less computational overhead. Therefore, (MOHAM) shows its validity and great potential for the
solution of nonlinear system of PDEs problems in science and engineering.

References

[1] M. Torrisi, R. Tracina, A. Valenti, A group analysis approach for a nonlinear differential
system arising in diffusion phenomena, J. Math. Phys. Vol. 37(9), (1996), 4758-4767.

[2] P.G. Drzain, R.S. Johnson, An introduction discussion the theory of solution and its diverse
applications, Camb. Uni. Pres. (1989).

[3] B. Abdel-Hamid, Exact solutions of some nonlinear evolution equations using symbolic
computations, Comp. Math. Appl. Vol. 40, (2000), 291-302.

[4] G. Bluman, S. Kumei, On the remarkable nonlinear diffusion equations, J. Math. Phys.,
Vol. 21(50), (1980), 1019-1023.

[5] LC. Chun, Fourier series based variational iteration method for a reliable treatment of
heat equations with variable coefficients, Int. J. Non-linr. Sci. Num. Simu., Vol. 10, (2007),
1383-1388.

[6] S.H. Chowdhury, A comparison between the Extended homotopy perturbation method and
adomian decomposition method for solving nonlinear heat transfer equations, J, Appl, Sci.,
Vol. 11(8), (2011), 1416-1420.



REFERENCES 551

[7] H. Yaghoobi, M. Tirabi, The application of differential transformation method to nonlinear
equation arising in heat transfer, Int. Comm. Heat Mass Transfer, Vol. 38(6), (2011), 815-
820.

[8] D.D. Ganji, The application of Hes homotopy perturbation method to nonlinear equation
arising in heat transfer, Phy. Lett. A, Vol. 355, (2006), 337-341.

[9] R. Bellman, Perturbation techniques in Mathematics, Phys. and Engg. Holt, Rinehart and
Winston, New York, (1964).

[10] J.D. Cole, Perturbation Methods in Applied Mathematics, Blaisedell, Waltham, MA. (1968).

[11] R.E. OMalley, Introduction to Singular Perturbation, Acad. Pres. New York. (1974).

[12] S.J. Liao, The proposed homotopy analysis technique for the solution of nonlinear problems,
PhD thesis, Shanghai Jiao Tong University, (1992).

[13] V. Marinca, N.Herisanu, I.Nemes, Optimal homotopy asymptotic method with application
to thin film flow, Central European Journal of Physics, Vol. 6, (2008), 64853.

[14] M. Sheikholeslami, Numerical investigation for CuO-H2O nanofluid flow in a porous chan-
nel with magnetic field using mesoscopic method, Journal of Molecular fluid, Vol. 249,
(20188), 73974.

[15] N. Herisanu, V, Marinca, Explicit analytical approximation to large amplitude nonlinear
osciallation of a uniform contilaver beam carrying an intermediate lumped mass and rotary
inertia, Meccanica, Vol. 45, (2010), 847-855.

[16] V. Marinca, N.Herisanu, An optimal homotopy asymptotic method applied to the steady flow
of a fourth-grade fluid past a porous plate, Applied Mathematics Letters, Vol. 22, (2009),
24551.

[17] V. Marinca, N.Herisanu, I.Nemes, A New Analytic Approach to Nonlinear Vibration of An
Electrical Machine, Proceeding of the Romanian Academy, Vol. 9, (2008) 229-236.

[18] V. Marinca, N.Herisanu, Determination of periodic solutions for the motion of a particle
on a rotating parabola by means of the Optimal Homotopy Asymptotic Method, Journal of
Sound and Vibration, Vol. 329(9), (2010), 1450-1459, doi: 10.1016/j.jsv.2009.11.005.

[19] H. Ullah, S. Islam, M. Fiza, An Extension of the Optimal Homotopy Asymptotic Method
to Coupled Schrdinger-KdV Equation, International journal of differential equations. Vol.
2014, Article ID 106934, 12 pages.

[20] H. Ullah, S. Islam, M. Idrees, M. Arif, Solution of Boundary Layer Problems with Heat
Transfer by Optimal Homotopy Asymptotic Method, Abstract and Applied Analysis Vol.
2013, Artcile ID 324869, 10 pages.

[21] H. Ullah, S. Islam, M. Idrees, M. Fiza, Application of Optimal Homotopy Asymptotic
Method to Doubly Wave Solutions of the Coupled Drinfeld-Sokolov-Wilson Equations, Math-
ematical Problem in Engineering Vol. 2013, Article ID 362816, 8 pages.

[22] H. Ullah, S. Islam, S. Sharidan, I. Khan, M. Fiza, Formulation and application of Optimal
Homotopy Asymptotic Method for coupled Differential Difference Equations, PLOS ONE
10.137/journal.pone.0120127, (2015).

[23] H. Ullah, S. Islam, I. Khan, L.C.C. Dennis, M. Fiza, Approximate Solution of Two-
Dimensional Nonlinear Wave Equation by Optimal Homotopy Asymptotic Method, Mathe-
matical Problems in Engineering Vol. 2015, Article ID 380104.



REFERENCES 552

[24] S. Abbasbandy, The application of homotopy analysis method to a generalized Hirota-
Satsuama coupled KdV equations, Phys. Lett. A., 361, (2007), 478-483.

[25] J. Liu, H. Li, Approximate analytic solutions of time fractional Hirota satsuam coupled KdV
and coupled MkdV equations, Abst. Appl. Anal. Vol. 2013, Arti. ID. 561980.