EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 17, No. 3, 2024, 1516-1538
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Numerical Finite-Difference Approximations of a
Coupled Reaction-Diffusion System with Gradient

Terms

Manar I. Khalil1,2,∗, Ishak Hashim1,3, Maan A. Rasheed4, Faieza Samat5,
Shaher Momani6,7

1 Department of Mathematical Sciences, Faculty of Science and Technology,
Universiti Kebangsaan Malaysiay, 43600 UKM Bangi Selangor, Malaysia
2 Department of Applied Geology, College of Science, Tikrit University, Iraq
3 Nonlinear Dynamics Research Center (NDRC), Ajman University, Ajman PO Box 346,
United Arab Emirates
4 Department of Mathematics,College of Basic Education, Mustansiriyah University, Bagh-
dad, Iraq
5 Pusat GENIUS@Pintar Negara, Universiti Kebangsaan Malaysiay, 43600 UKM Bangi
Selangor, Malaysia
6 Nonlinear Dynamics Research Center (NDRC), Ajman University, Ajman PO Box 346,
United Arab Emirates
7 Department of Mathematics,Faculty of Science, Jordan University, Amman 11942, Jordan

Abstract. This study focuses on the derivation of explicit and implicit finite difference formu-
las.The objective of this study is to derive an estimation of the blow-up time for a coupled reaction-
diffusion system incorporating gradient terms, employing numerical finite difference approxima-
tions. Furthermore, an examination is conducted on the consistency, stability, and convergence
of the proposed schemes. Additionally, the study presents two numerical experiments. In each
instance, the numerical blow-up time is calculated benefit the suggested methodologies, employing
varying space steps and non-fixed time-stepping. The numerical findings obtained demonstrate
that the blow-up time sequence exhibits convergence as the space step decreases. Moreover, the
numerical orders of convergence for the blow-up time goes well with the theoretical orders observed
in the numerical solutions.

2020 Mathematics Subject Classifications: 47H05, 90C33

Key Words and Phrases: Explicit and Implicit finite difference formulas, blow-up time, Con-
vergence analysis, consistency, stability, gradient terms

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v17i3.5246

Email addresses: p119266@siswa.ukm.edu.my (Manar I. Khalil),
xxxxx@email.com (MISSING AUTHOR EMAILS)

https://www.ejpam.com 1516 © 2024 EJPAM All rights reserved.



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1517

1. Introduction

There are many non-linear partial differential equations of the parabolic type that
cannot be extended globally in time for a given initial data and become unbounded in
a limited time. The occurrence of this phenomenon is referred to as blow-up, and it
can manifest in time-dependent nonlinear equations, provided that the non-linear factors
possess sufficient strength, as evidenced by references [9, 10, 14, 17, 21]. This phenomenon
can be expanded to coupled systems when each dependent variable has a finite increase
within a specific time-frame. In this particular scenario, it is asserted that the dependent
variables exhibit synchronous increases [7, 18, 19, 22, 23].

In this paper, we consider a coupled parabolic system:

ut = uxx − |ux|q1 + |v|p1 , vt = vxx − |vx|q2 + |u|p2 , (x, t) ∈ (0, 1)× (0, T ),
u(0, t) = u(1, t) = 0,
v(0, t) = v(1, t) = 0, t ∈ (0, T ),
u(x, 0) = u0(x), v(x, 0) = v0(x), x ∈ (0, 1),

 (1.1)

where p1, p2 > 1, q1, q2 ∈ (1, 2), u0, v0 ≥ 0. For some details regarding the local existence
and uniqueness and blow-up of this system, see [21]. The system (1.1) has been studied
in [15] ,where Ω ∈ R3 and is a bounded convex domain, and

p = p1 = p2,

(
α = β =

1

1− p

)
,

q = q1 = q2; p > q > 1.

This is shown as the traditional solution of the system blows up (becomes unbounded) in
the W-norm, where

W (t) =

∫
Ω

(
u2p + v2p

)
dx,

then blow-up time for this problem can be estimated from below as follows:

T ≥ 1

2AW 2
0

,

where W (0) = W0 =
∫
Ω

(
u2p0 + v2p0

)
dx and A is a constant which depends on the data.

In [21], with some restricted conditions on system (1.1), it was shown that the upper
blow-up rate estimates for this solution and its gradients terms take the following forms:

u(x, t) ≤ C1(T − t)−α, |∇u(x, t)| ≤ C1(T − t)−
(1+2α)

2 ,

v(x, t) ≤ C2(T − t)−β, |∇v(x, t)| ≤ C2(T − t)−
(1+2β)

2 ,

where (x, t) ∈ Ω× (0, T ), C1, C2 > 0, α = p1+1
p1p2−1 and β = p2+1

p1p2−1 .
The numerical approximation of time-dependent parabolic problems has been studied

by many authors, see for instance [1–4, 11–13, 16, 23, 24]. The main aim of these works



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1518

was to estimate the numerical blow-up time of the considered problems. Blow-up means
the solution of the problem becomes Unbounded in finite time. The semi-linear coupled
reaction-diffusion system was one of the investigated problems :

ut = uxx + vp, vt = vxx + uq, x ∈ (0, 1), t ∈ (0, T ),
u(0, t) = u(1, t) = 0, v(0, t) = v(1, t) = 0, t ∈ (0, T ),
u(x, 0) = u0(x), v(x, 0) = v0(x), x ∈ (0, 1),

 (1.2)

where p, q > 1 and u0(x), v0(x) are both smooth functions.
The semi-discrete approximation problem for system (1.2) is derived in [13], and the-

oretical notions regarding the rapprochement of blow-up times and blow-up solutions of
the semi-discrete issue are proven . Furthermore, this study introduces two alternative ap-
proximations to problem (1.2), namely Euler schemes, which incorporate a non-fixed time
stepping formula. Furthermore, two numerical experiments were conducted to provide em-
pirical support for the numerical results. Specifically, the authors conducted estimations
for a numerical blow-up time, error bounds, CPUt, and numerical order of convergence.
The Chipot-Weissler equation [5]: refers to the single equation of system (1.1).

ut = uxx + up − |ux|q , (x, t) ∈ (0, 1)× (0, T ),
u(0, t) = u(1, t) = 0, t ∈ (0, T ),
u(x, 0) = u0(x), x ∈ (0, 1).

 (1.3)

System (1.3) is associated with a widely recognized model that has emerged in the field of
population dynamics research [20]. It is worth noting that the inclusion of the term may be
gradation offer a damping effect that mitigates blow-up. Consequently, numerous scholars
have expressed interest in investigating the impact of the gradient term on blow-up char-
acteristics, including blow-up set and blow-up rate estimations [6, 8, 25]. The problem’s
numerical approximation is examined in reference [24]. The semi-discrete approximation
of equation (1.1) was derived by the authors, who also conducted theorems to establish the
convergence and blow-up of the semi-discrete problem. Moreover, the authors put forth
two approximation equations for (1.1) that are fully discrete: explicit Euler equations and
implicit finite difference equations. These equations are formulated using a non-fixed time-
stepping algorithm. Furthermore, two numerical experiments are conducted to determine
the numerical blow-up time, error limits, and numerical convergence order. The primary
objective of this study is to employ numerical finite difference approximations in order to
estimate the blow-up time of a coupled reaction-diffusion system that incorporates gra-
dient terms. The explicit and implicit finite difference formulas are developed in parts
two and three, respectively. Furthermore, an examination is conducted on the stability
,convergence , and consistency of the proposed schemes. Two numerical experiments are
presented in a fourth part. In each instance, the numerical blow-up time is calculated
utilizing the suggested methodologies, employing varying space stages and non-fixed time
steps. In the final section, we present a set of conclusions and potential avenues for future
works.



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1519

2. Explicit Euler scheme

Here, we derive the explicit fully-discrete finite difference formulas for the problem
(1.1),by approximating the time derivative in problem(1.2), using the forward finite dif-
ference formula:
We consider that Un

i and V n
i , are the approximate values of u (xi, tn) and v (xi, tn), re-

spectively, where tn+1 = tn + kn, I be a positive integer, and consider the grid: xi = ih ,
0 ≤ i ≤ I, h = 1

I , tn+1 = tn + kn, xi+1 = xi + h, h is the space-step,and kn is the time-step
.
We approximate ut, vt by forward finite difference formulas at the mesh-point (xi, tn), as
below :

ut|ni =
Un+1
i − Un

i

kn
+O (kn) , vt|ni =

V n+1
i − V n

i

kn
+O (kn) ,

while uxx, vxx The second-order central finite difference formula is used to approximate
the following:

uxx|ni =
Un+1
i − 2Un

i + Un
i−1

h2
+O

(
h2
)
, vxx|ni =

V n+1
i − 2V n

i + V n
i−1

h2
+O

(
h2
)
.

Also, ux, vx are approximated using first-order central finite difference formula as follows:

ux|ni =
Un
i+1 − Un

i−1

2h
+O

(
h2
)
, vx|ni =

V n
i+1 − V n

i−1

2h
+O

(
h2
)
.

Substituting in system (1.1) yields

Un+1
i − Un

i

kn
=

Un
i − 2Un

i + Un
i−1

h2
+ F (δxU

n
i , V

n
i ) ,

V n+1
i − V n

i

kn
=

V n
i − 2V n

i + V n
i−1

h2
+G (δxV

n
i , Un

i ) ,

where

F (δxU
n
i , V

n
i ) = (V n

i )p1 −
∣∣∣∣Un

i+1 − Un
i−1

2h

∣∣∣∣q1 , G (δxV
n
i , Un

i ) = (Un
i )

p2 −
∣∣∣∣V n

i+1 − V n
i−1

2h

∣∣∣∣q2 .
These equations can be rewritten as

Un+1
i = (1− 2rn)U

n
i + rn

(
Un
i+1 + Un

i−1

)
+ knF (δxU

n
i , V

n
i ) , (2.1)

V n+1
i = (1− 2rn)V

n
i + rn

(
V n
i+1 + V n

i−1

)
+ knG (δxV

n
i , Un

i ) , (2.2)

where

i = 1, 2, . . . , I − 1, n = 0, 1, 2, . . . ,

Un
h =

(
Un
1 , U

n
2 , . . . , U

n
1−I

)
, V n

h =
(
V n
1 , V n

2 , . . . , V n
1−I

)
, rn =

kn
h2

.



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1520

Furthermore, to ensure that the stability (convergence) criterion of the explicit scheme is
satisfied, the non-stationary time step formula is taken as below :

kn = min

(
h2

3
,

hα∥∥Un
h

∥∥ , hα∥∥Vn
h

∥∥
)
, α ≥ 1. (2.3)

Lemma 1. The functions F,G satisfy Lipshtiz condition, that is there exists positive
constants L1, L2, L3, L4 > 0, such as∣∣∣F (δxU

n
i , V

n
i )− F

(
δxŨ

n
i , Ṽ

n
i

)∣∣∣ ≤ L1

∣∣∣δx (Un
i − Ũn

i

)∣∣∣+ L2

∣∣∣V n
i − Ṽ n

i

∣∣∣ , (2.4)∣∣∣G (Un
i , δxV

n
i )−G

(
Ũn
i , δxṼ

n
i

)∣∣∣ ≤ L3

∣∣∣δx (V n
i − Ṽ n

i

)∣∣∣+ L4

∣∣∣Un
i − Ũn

i

∣∣∣ , (2.5)

where Un
i , V

n
i , Ũn

i , Ṽ
n
i (i = 0, 1, 2, . . . I) are bounded values.

Proof. We denote f(s) = sp and we use the mean value theorem we get∣∣∣(V n
i )p1 −

(
Ṽ n
i

)p1∣∣∣ ≤ ∣∣f ′ (Zi)
∣∣ ∣∣∣V n

i − Ṽ n
i

∣∣∣ , (2.6)

where Zi is an intermediate value between V n
i and V̄ n

i . Since V n
h , V n

h are bounded, then
there exists C1 > 0 such that |f ′ (Zi)| < C1. Now, we denote g(s) = |s|q1 , with using the
mean value theorem, we find∣∣∣∣∣
∣∣∣∣Un

i+1 − Un
i−1

2h

∣∣∣∣q1 −
∣∣∣∣∣ Ũn

i+1 − Ũn
i−1

2h

∣∣∣∣∣
q1
∣∣∣∣∣ ≤ ∣∣g′ (εi)∣∣

∣∣∣∣∣
(
Un
i+1 − Un

i−1

2h

)
−

(
Ũn
i+1 − Ũn

i−1

2h

)∣∣∣∣∣ , (2.7)
where εi is an intermediate value between

a1i =

(
Un
i+1 − Un

i−1

2h

)
and a2i =

(
Ũn
i+1 − Ũn

i−1

2h

)
,

since Un
h Ũ

n
h are bounded in [0, T ]. Therefore, there exists C2 > 0 such that |g′ (εi)| ≤ C2.

From (2.6) and (2.7), we get (2.4). Similarly, we can show that (2.5) is held.

2.1. The algorithm steps for Euler Explicit method

1. Input h, U0
h , V

0
h ,p1,p2, q1, q2, α

2. Put n = 0;

3. Choose kn according to (2.3).

4. Compute the numerical vectors: Un+1
h , V n+1

h , using the explicit formula (2.1) and
(2.2).

5. For n = 1, 2, . . . .. , repeat steps 3,4 until for n = m, we get ∥Un
h ∥∞ ≥ 1015, or

∥V n
h ∥∞ ≥ 1015

6. The numerical blow-up time is tm =
∑m

n=0 kn.



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1521

2.2. Local Truncation error of Explicit Euler scheme

Theorem 1. Let (Tn
ui, T

n
vi) be The local error truncation of the explicit Euler chart (2.1)

and (2.2) at the grid-point (xi, tn). after that, Here it isC1, C2, C3, C4 > 0 like this

|Tn
ui| ≤ C1k + C2h

2, |Tn
vi| ≤ C3k + C4h

2,

where k = maxn∈N kn, i.e. Tn
ui = O

(
k + h2

)
and Tn

vi = O
(
k + h2

)
.

Proof. Replace the precise solution uni = u (xi, tn) and vni = v (xi, tn) into the explicit
(2.1) yields

Tn
ui =

(
un+1
i − uni

)
− kn

h2
(
uni+1 − 2uni + uni−1

)
− kn (v

n
i )

p1 + kn

∣∣∣∣Un
i+1 − Un

i−1

2h

∣∣∣∣q1
= kn

[
∂uni
∂t

+O(kn)

]
− kn

[
∂2uni
∂x2

+O
(
h2
)]

− kn (v
n
i )

p1 + kn

[∣∣∣∣∂uni∂x

∣∣∣∣q1 +O
(
h2
)]

= kn

[
∂uni
∂t

− ∂2uni
∂x2

− (vni )
p1 +

∣∣∣∣∂uni∂x

∣∣∣∣q1]+ kn
[
O(kn) +O

(
h2
)]

.

From system (1.1), assuming all the partial derivatives are bounded at the grid-point
(xi, tn), we obtain |Tn

ui| ≤ C1kn + C2h
2 ≤ C1k + C2h

2, where C1, C2 > 0, i.e. Tn
ui =

O
(
k + h2

)
, C > 0.

Likewise, we show that Tn
vi = O

(
k + h2

)
, C > 0.

2.3. The stability analysis of explicit Euler method

The stability condition for explicit Euler scheme can be obtained based on the the
following definition

Definition 1. [13]. For i = 1, 2, . . . , I − 1, set En
u = (enui), E

n
v = (envi), e

n
ui = uni − Un

i

and envi = vni − V n
i , where uni = u(xi, tn), v

n
i = v(xi, tn) and (Un

i , V
n
i ) are a accurate and

numerical solutions of a one-dimensional time-dependent coupled system of two PDEs,
Consecutively. For any hypothetical initial rounding error (E0

u, E
0
v), we say that the nu-

merical solution is stable, If a positive number is foundµ freelance on the space-step (h)
and time-step (k), as follows

∥En
u∥ ≤ µmax

{∥∥E0
u

∥∥ ,∥∥E0
v

∥∥} and ∥En
v ∥ ≤ µmax

{∥∥E0
u

∥∥ ,∥∥E0
v

∥∥} ,
where ∥En

u∥ = max1≤i≤I |enui| and ∥En
v ∥ = max1≤i≤I |envi| , n = 0, 1, 2, . . . .

Theorem 2. The explicit formulae(2.1), (2.2) is stable, if (1 − 2r) ≥ 0, where k =
maxn∈N kn, r = k/h2.



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1522

Proof. For proof this theorem, we use the maximum error stability-technicality [27].
Let enui = uni − Un

i , envi = vni − V n
i , (uni = u(xi, tn), v

n
i = v(xi, tn)) be Accurate

solution to the issue (1.1).
To finish this, we implement mathematical induction. For n = 1, set

∥∥E1
u

∥∥ = max1≤i≤I

∣∣e1ui∣∣ =∣∣∣e1uj∣∣∣ , ∥∥E1
v

∥∥ = max1≤i≤I

∣∣e1vi∣∣ = ∣∣e1um∣∣. By substituting enuj in the explicit formula (2.1),

we obtain

e1uj = (1− 2r1)e
0
uj + r1

(
e0uj+1 + e0uj−1

)
+ k1

(
F
(
δxu

0
j , v

0
j

)
− F

(
δxU

0
j , V

0
j

))
.

Since (1− 2r1) ≥ 0, we have∣∣e1uj∣∣ ≤ (1− 2r1)
∣∣e0uj∣∣+ r1

(∣∣e0uj+1

∣∣+ ∣∣e0uj−1

∣∣)+ k1
∣∣F (δxu0j , v0j )− F

(
δxU

0
j , V

0
j

)∣∣ .
From (1.2), it follows that∣∣e1uj∣∣ ≤ (1− 2r1)

∣∣e0uj∣∣+ r1
(∣∣e0uj+1

∣∣+ ∣∣e0uj−1

∣∣)+ k1L1

∣∣δx (u0j − U0
j

)∣∣+ k1L2

∣∣v0j − V 0
j

∣∣∣∣e1uj∣∣ ≤ (1− 2r1)
∥∥E0

u

∥∥+ 2r1
∥∥E0

u

∥∥+ k1L1

∣∣δxe0uj∣∣+ k1L2

∣∣e0vj∣∣∥∥E1
u

∥∥ ≤
∥∥E0

u

∥∥+ 2k1L1

∥∥E0
u

∥∥+ k1L2

∥∥E0
v

∥∥
≤ (1 + 2k1L1 + k1L2)Max

{∥∥E0
u

∥∥ , ∥∥E0
v

∥∥}
≤ (1 + 3kL)Max

{∥∥E0
u

∥∥ ,∥∥E0
v

∥∥} .
Similarly, we can show that∥∥E1

v

∥∥ ≤ (1 + 2kL3 + kL4)Max
{∥∥E0

u

∥∥ , ∥∥E0
v

∥∥}
≤ (1 + 3kL)Max

{∥∥E0
u

∥∥ , ∥∥E0
v

∥∥} ,
where L = Max {L1, L2, L3, L4}. Now, we suppose that

∥Es
u∥ ≤ (1 + 3kL)sMax

{∥∥E0
u

∥∥ ,∥∥E0
v

∥∥} , s = 1, 2, 3, . . . n,

∥Es
v∥ ≤ (1 + 3kL)sMax

{∥∥E0
u

∥∥ ,∥∥E0
v

∥∥} , s = 1, 2, 3, . . . n.

For n + 1, set
∥∥En+1

u

∥∥ = max1≤i≤I

∣∣en+1
ui

∣∣ = ∣∣∣en+1
uj

∣∣∣ ,∥∥En+1
v

∥∥ = max1≤i≤l

∣∣en+1
vi

∣∣ = ∣∣en+1
um

∣∣.
By substituting enuj in the explicit formula (2.1), we obtain

en+1
uj = (1− 2rn)e

n
uj + rn

(
enuj+1 + enuj−1

)
+ kn

(
F
(
δxu

n
j , v

n
j

)
− F

(
δxU

n
j , V

n
j

))
.

Since (1− 2r) ≥ 0, we have∣∣∣en+1
uj

∣∣∣ ≤ (1− 2rn)
∣∣enuj∣∣+ rn

(∣∣enuj+1

∣∣+ ∣∣enuj−1

∣∣)+ kn
∣∣F (δxunj , vnj )− F

(
δxU

n
j , V

n
j

)∣∣ .
Thus, from (1.2), it follows that∣∣∣en+1

uj

∣∣∣ ≤ (1− 2rn) ∥En
u∥+ 2rn ∥En

u∥+ knL1

∣∣δx (unj − Un
j

)∣∣+ knL2

∣∣vnj − V n
j

∣∣



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1523

≤ ∥En
u∥+ knL1

∣∣δxenuj∣∣+ knL2

∣∣envj∣∣∥∥En+1
u

∥∥ ≤ ∥En
u∥+ 2knL1 ∥En

u∥+ knL2 ∥En
v ∥

≤ (1 + 3kL)Max {∥En
u∥ , ∥En

v ∥}
≤ (1 + 3kL)(1 + 3kL)nMax

{∥∥E0
u

∥∥ ,∥∥E0
v

∥∥}
≤ (1 + 3kL)n+1Max

{∥∥E0
u

∥∥ , ∥∥E0
v

∥∥}
≤ exp(3(n+ 1)kL)Max

{∥∥E0
u

∥∥ ,∥∥E0
v

∥∥}
= exp (3tn+1L)Max

{∥∥E0
u

∥∥ ,∥∥E0
v

∥∥} .
Likewise, ∥∥En+1

v

∥∥ ≤ exp (3tn+1L)Max
{∥∥E0

u

∥∥ , ∥∥E0
v

∥∥} .
So, this implies the stability for explicit scheme, if r ≤ 1/2.

2.4. Convergence analysis of explicit Euler scheme

Theorem 3. The explicit Euler formulas (2.1), (2.2) are convergent with: O
(
k + h2

)
, if

(1− 2r) ≥ 0, where k = maxn∈N kn, r = k/h2.

Proof. Set enui = uni −Un
i , e

n
vi = vni −V n

i , (uni = u(xi, tn), v
n
i = v(xi, tn)) is the accurate

solution to the problem (1.1). suppose that e0ui = 0, e0vi = 0, ∀i = 0, 1, . . . , I. We seek for
C > 0, such that

en+1
ui ≤ C

(
k + h2

)
, en+1

vi ≤ C
(
k + h2

)
, n = 0, 1, . . . .

To finish this, We use the Tech of mathematical induction For n = 1, we set
∣∣∣e1uj∣∣∣ =

max1≤i≤I−1

∣∣e1ui∣∣. Substituting e1uj in the explicit formula (2.1) yields that

e1uj = (1− 2rn)e
0
uj + rn

(
e0uj+1 + e0uj−1

)
+ kn

(
F
(
δxu

0
j , v

0
j

)
− F

(
δxU

0
j , V

0
j

))
+ Tui

0

= kn
(
F
(
δxu

0
j , v

0
j

)
− F

(
δxU

0
j , V

0
j

))
+ Tui

0.

From Lemma (1) we obtain∣∣e1uj∣∣ ≤ knL1

∣∣δx (u0j − U0
j

)∣∣+ knL2

∣∣v0j − V 0
j

∣∣+ ∣∣T 0
ui

∣∣ = ∣∣T 0
ui

∣∣ ≤ C
(
k + h2

)
.

Hence
∣∣e1ui∣∣ ≤ C

(
k + h2

)
, i = 1, 2, . . . I − 1. Likewise, we show that∣∣e1vi∣∣ ≤ C

(
k + h2

)
, i = 1, 2, . . . , I − 1.

assume
∣∣eSui∣∣ ≤ Cs

(
k + h2

)
,
∣∣eSvi∣∣ ≤ Cs

(
k + h2

)
, s = 0, 1, 2, . . . , n, Cs > 0. Let C∗ =

max0≤s≤nCs. For n+ 1, we set∣∣∣en+1
uj

∣∣∣ = max1≤i≤I−1

∣∣en+1
ui

∣∣.
Replace en+1

uj in the explicit scheme (2.1) yields

en+1
uj = (1− 2rn)e

n
uj + rn

(
enuj+1 + enuj−1

)
+ kn

(
F
(
δxu

n
j , v

n
j

)
− F

(
δxU

n
j , V

n
j

))
+ Tn

ui.



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1524

Thus∣∣∣en+1
uj

∣∣∣ ≤ (1− 2rn) ∥En
u∥+ 2rn ∥En

u∥+ kn
∣∣F (δxunj , vnj )− F

(
δxU

n
j , V

n
j

)∣∣+ ∣∣Tn
uj

∣∣
≤ (1− 2rn) ∥En

u∥+ 2rn ∥En
u∥+ knL1

∣∣δx (unj − Un
j

)∣∣+ knL2

∣∣vnj − V n
j

∣∣+ ∣∣Tn
uj

∣∣ .
From Lemma (1), we have∥∥En+1

u

∥∥ ≤ ∥En
u∥+ knL1

∣∣δxenuj∣∣+ knL2

∣∣envj∣∣+ ∣∣Tn
j

∣∣
≤ ∥En

u∥+ 2knL1 ∥En
u∥+ knL2 ∥En

v ∥+
∣∣Tn

j

∣∣
≤ (1 + 2knL1 + knL2)C

∗ (k + h2
)
+ C

(
k + h2

)
= [(1 + 2kL1 + kL2)C

∗ + C]
(
k + h2

)
.

Next
∥∥En+1

u

∥∥ ≤ C
(
k + h2

)
, n = 0, 1, . . .. Likewise, we show that∥∥En+1
v

∥∥ ≤ C
(
k + h2

)
, n = 0, 1, . . . .

3. Implicit Euler scheme

Now, in order to obtain the implicit fully-discrete finite difference formulae for problem
(1.1), we will employ the backward finite difference formula to approximate the time
derivative in problem (1.1). In this analysis, we define Un

i and V n
i , as the estimated values

of u (xi, tn) and v (xi, tn), consecutively, let tn+1 = tn + kn, where I represent a positive
integer, and proceed to examine the grid:: xi = ih, 0 ≤ i ≤ I, h = 1

I , tn+1 = tn + kn,
xi+1 = xi + h, h is the space-step, kn is the time-step.

At the grid point (xi, tn+1),we employ backward finite difference formulas to approxi-
mate the values of ut, vt as below :

∂u

∂t

∣∣∣∣n+1

i

=
1

kn

(
Un+1
i − Un

i

)
+O (kn) ,

∂v

∂t

∣∣∣∣n+1

i

=
1

kn

(
V n+1
i − V n

i

)
+O (kn) .

While uxx, vxx are approximated as using the second-order central finite difference formulas
below:

uxx|n+1
i =

Un+1
i − 2Un+1

i + Un+1
i−1

h2
+O

(
h2
)
,

vxx|n+1
i =

V n+1
i − 2V n+1

i + V n+1
i−1

h2
+O

(
h2
)
.

Moreover, the non-linear terms are approximated as follows:

F
(
δxU

n+1
j , V n+1

j

)
= F

(
δxU

n
j , V

n
j

)
+O(k),

G
(
δxV

n+1
i , Un+1

i

)
= G (δxV

n
i , Un

i ) +O(k).



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1525

Substituting all of these formulas into system (1.1) gives

Un+1
i − Un

i

kn
=

Un+1
i+1 − 2Un+1

i + Un+1
i−1

h2
+ F

(
δxU

n
j , V

n
j

)
, 1 ≤ i ≤ I − 1,

V n+1
i − V n

i

kn
=

V n+1
i+1 − 2V n+1

i + V n+1
i−1

h2
+G (δxV

n
i , Un

i ) , 1 ≤ i ≤ I − 1.

Thus

(1 + 2rn)U
n+1
i − rn

(
Un+1
i+1 + Un+1

i−1

)
= Un

i + knF
(
δxU

n
j , V

n
j

)
, (3.1)

(1 + 2rn)V
n+1
i − rn

(
V n+1
i+1 + V n+1

i−1

)
= V n

i + knG (δxV
n
i , Un

i ) , (3.2)

where rn = kn
h2 , 1 ≤ i ≤ I − 1, n = 0, 1, 2, . . ..

Moreover, The non-fixed time-stepping formula is employed due to the lack of condition
stability in the implicit Euler chart.

kn = min

(
hα∥∥Un
h

∥∥ , hα∥∥Vn
h

∥∥
)
, α ≥ 1. (3.3)

We write (3.1)–(3.2) in matrix form as below:

(I − rnhH)Un+1
h = Un

h + knF (V n
h ) , (3.4)

(I − rnhH)V n+1
h = V n

h + knG (Un
h ) , (3.5)

where

H =


−2 1 0 · · · 0
1 −2 1 · · · 0

. . .
...

. . .

0 · · · 0 1 −2


(I−1)×(I−1)

,

F (V n
h ) =

(
f (V n

1 ) , f (V n
2 ) , . . . , f

(
V n
I−1

))T
, G (Un

h ) =
(
g (Un

1 ) , g (U
n
2 ) , . . . , g

(
Un
I−1

))T
.

3.1. The algorithm steps for Euler implicit method

1. Input h, U0
h , V

0
h ,p1,p2, q1, q2, α

2. Put n = 0;

3. Choose kn according to (3.3).

4. Compute the numerical vectors: Un+1
h , V n+1

h , by solving the linear systems (3.4) and
(3.5).

5. For n = 1, 2, . . . .. , repeat steps 3,4 until for n = m, we get ∥Un
h ∥∞ ≥ 1015, or

∥V n
h ∥∞ ≥ 1015



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1526

6. The numerical blow-up time is tm =
∑m

n=0 kn.

7.

3.2. Local Truncation error of implicit Euler scheme

Theorem 4. Let (Tn
ui, T

n
vi) be the local truncation error of the implicit Euler scheme (3.1)–

(3.2) at the mesh point (xi, tn+1). Then, there exist C1, C2, C3, C4 > 0, such that
∣∣Tn+1

ui

∣∣ ≤
C1k+C2h

2 and
∣∣Tn+1

vi

∣∣ ≤ C3k+C4h
2, where k = maxn∈N kn, i.e. Tn+1

ui = O
(
k + h2

)
and

Tn+1
vi = O

(
k + h2

)
.

Proof. Replace the precise solution uni = u (xi, tn) , v
n
i = v (xi, tn) into the implicit

Euler scheme (3.1) yields

Tn+1
ui =

(
un+1
i − uni

)
− kn

h2
[
un+1
i+1 − 2un+1

i + un+1
i−1

]
− knF

(
δxU

n
j , V

n
j

)
= kn

[
∂un+1

i

∂t
+O(kn)

]
− kn

[
∂2un+1

i

∂x2
+O

(
h2
)]

= −k

[(
vn+1
i

)p1 − ∣∣∣∣∂un+1
i

∂x

∣∣∣∣q1 +O(kn) +O
(
h2
)]

.

It follows that

Tn+1
ui = kn

[
∂un+1

i

∂t
−

∂2un+1
i

∂x2
−
(
vn+1
i

)p1
+

∣∣∣∣∂un+1
i

∂x

∣∣∣∣q1
]
+ kn

[
O(kn) +O

(
h2
)]

.

From equation (1.1) and Presuming all the partial derivatives are bounded at the mesh-
point (xi, tn+1), we obtain

∣∣Tn+1
ui

∣∣ ≤ C1kn + C2h
2 ≤ C1k + C2h

2, where C1, C2 > 0. , i.e.∣∣Tn+1
i

∣∣ = O
(
k + h2

)
, C > 0. Similarly, we show that there exist C3, C4 > 0 such that∣∣Tn+1

vi

∣∣ ≤ C3k + C4h
2.

3.3. The stability analysis of implicit Euler method

Theorem 5. The implicit Euler scheme (3.1)–(3.2) is unconditionally stable.

Proof. To demonstrate the validity of this theorem, the maximum error stability
technique can be employed. [6].

Let enui = uni − Un
i , envi = vni − V n

i , (uni = u (xi, tn) , v
n
i = v (xi, tn)) is the accurate

solution to the problem (1.1).
To finish this, We implement the mathematical induction. For n = 1 and setting∥∥E1
u

∥∥ = max1≤i≤I

∣∣e1ui∣∣ = ∣∣∣e1uj∣∣∣, k = maxn∈N kn, r = k/h2, we have∣∣e1uj∣∣ = (1 + 2r1)
∣∣e1uj∣∣− r1

(∣∣e1uj∣∣+ ∣∣e1uj∣∣)



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1527

≤ (1 + 2r1)
∣∣e1uj∣∣− r1

(∣∣e1uj+1

∣∣+ ∣∣e1uj−1

∣∣)
≤

∣∣(1 + 2r1)e
1
uj − r1

(
e1uj+1 + e1uj−1

)∣∣
=

∣∣e0uj + k1
(
F
(
δxu

0
j , v

0
j

)
− F

(
δxU

0
j , V

0
j

))∣∣ .
Lemma (1) then gives∣∣e1uj∣∣ ≤ ∣∣e0j ∣∣+ k1

∣∣F (δxu0j , v0j )− F
(
δxU

0
j , V

0
j

)∣∣
≤
∣∣e0j ∣∣+ k1L1

∣∣δx (u0j − U0
j

)∣∣+ k1L2

∣∣v0j − V 0
j

∣∣
=
∣∣e0j ∣∣+ k1L1

∣∣δxe0uj∣∣+ k1L2

∣∣e0vj∣∣ .
It follows that∥∥E1

u

∥∥ ≤
∥∥E0

u

∥∥+ 2kL1

∥∥E0
u

∥∥+ kL2

∥∥E0
v

∥∥ ≤ (1 + 3kL)Max
{∥∥E0

u

∥∥ , ∥∥E0
v

∥∥} .
Similarly, we can show that∥∥E1

v

∥∥ ≤
∥∥E0

v

∥∥+ 2kL3

∥∥E0
v

∥∥+ kL4

∥∥E0
u

∥∥ ≤ (1 + 3kL)Max
{∥∥E0

u

∥∥ , ∥∥E0
v

∥∥} ,
where L = Max {L1, L2, L3, L4}.

Now, we suppose that

∥Es
u∥ ≤ (1 + 3kL)sMax

{∥∥E0
u

∥∥ , ∥∥E0
v

∥∥} , s = 1, 2, 3, . . . n,

∥Es
v∥ ≤ (1 + 3kL)sMax

{∥∥E0
u

∥∥ , ∥∥E0
v

∥∥} , s = 1, 2, 3, . . . n.

For n+1 and setting
∥∥En+1

u

∥∥ = max1≤i≤I

∣∣en+1
ui

∣∣ = ∣∣∣en+1
uj

∣∣∣ and ∥∥En+1
v

∥∥ = max1≤i≤I

∣∣en+1
vi

∣∣,
k = maxn∈N kn, r = k/h2 yields∣∣∣en+1

uj

∣∣∣ = (1 + 2rn)
∣∣∣en+1

uj

∣∣∣− rn

(∣∣∣en+1
uj

∣∣∣+ ∣∣∣en+1
uj

∣∣∣)
≤ (1 + 2rn)

∣∣∣en+1
uj

∣∣∣− rn

(∣∣∣en+1
uj+1

∣∣∣+ ∣∣∣en+1
uj−1

∣∣∣)
≤
∣∣∣(1 + 2rn)e

n+1
uj − rn

(
en+1
uj+1 + en+1

uj−1

)∣∣∣
=
∣∣enuj + kn

(
F
(
δxu

n
j , v

n
j

)
− F

(
δxU

n
j , V

n
j

))∣∣ .
By Lemma (1), we obtain∣∣∣en+1

uj

∣∣∣ ≤ ∣∣enuj∣∣+ kn
∣∣fF (δxunj , vnj )− F

(
δxU

n
j , V

n
j

)∣∣
≤
∣∣enuj∣∣+ knL1

∣∣δx (unj − Un
j

)∣∣+ knL2

∣∣vnj − V n
j

∣∣
=
∣∣enj ∣∣+ knL1

∣∣δxenuj∣∣+ knL2

∣∣envj∣∣ .
It follows that∥∥En+1

u

∥∥ ≤ ∥En
u∥+ 2knL ∥En

u∥+ knL ∥En
v ∥



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1528

≤ (1 + 3kL)nMax
{∥∥E0

u

∥∥ , ∥∥E0
v

∥∥}+ 3kL(1 + 3kL)nMax
{∥∥E0

u

∥∥ ,∥∥E0
v

∥∥}
≤ (1 + 3kL)n+1Max

{∥∥E0
u

∥∥ ,∥∥E0
v

∥∥}
≤ exp(3(n+ 1)kL)Max

{∥∥E0
u

∥∥ , ∥∥E0
v

∥∥}
= exp (3tn+1L)Max

{∥∥E0
u

∥∥ ,∥∥E0
v

∥∥} .
Similarly, we can show that∥∥En+1

v

∥∥ ≤ exp (3tn+1L)Max
{∥∥E0

u

∥∥ , ∥∥E0
v

∥∥} . where L = Max {L1, L2, L3, L4}

So, the implicit Euler scheme (2.4)–(2.5) is unconditionally stable.

3.4. Convergence analysis of implicit Euler scheme

Theorem 6. The implicit Euler formula (3.1)–(3.2) is convergent with r = k/h2 > 0, if
r > 0, where k = maxn∈N kn.

Proof. Set enui = uni − Un
i , e

n
vi = vni − V n

i , and let (uni = u(xi, tn) and vni = v(xi, tn))
be Accurate solution to the issue (1.1). suppose that e0ui = 0, e0vi = 0, ∀i = 0, 1, . . . , I. We
aim to show that there exists C > 0 such that

en+1
ui ≤ C

(
k + h2

)
, en+1

vi ≤ C
(
k + h2

)
, n = 0, 1, . . . .

To finish this, we use the mathematical induction tech. For n = 1, we set |e1uj | =

max1≤i≤I−1 |e1ui|. Thus

|e1uj | = (1 + 2r1)|e1uj | − r1(|e1uj |+ |e1uj |)
≤ (1 + 2r1)

∣∣e1uj∣∣− r1
(∣∣e1uj+1

∣∣+ ∣∣e1uj−1

∣∣)
≤
∣∣(1 + 2r1)e

1
uj − r1

(
e1uj+1 + e1uj−1

)∣∣
=
∣∣e0uj + k1

(
F
(
δxu

0
j , v

0
j

)
− F

(
δxU

0
j , V

0
j

))
+ T 0

i

∣∣
=
∣∣k1 (F (δxu0j , v0j )− F

(
δxU

0
j , V

0
j

))
+ T 0

j

∣∣ .
From Lemma (1), we have∣∣e1uj∣∣ ≤ k1L1

∣∣δx (u0j − U0
j

)∣∣+ k1L2

∣∣v0j − V 0
j

∣∣+ ∣∣T 0
j

∣∣ ,
= k1L1

∣∣δxe0uj∣∣+ k1L2

∣∣e0vj∣∣+ ∣∣T 0
j

∣∣
where L = Max {L1, L2}. So,∣∣e1uj∣∣ ≤ 2kL1

∥∥E0
u

∥∥+ kL2

∥∥E0
v

∥∥+ ∣∣T 0
i

∣∣ = ∣∣T 0
i

∣∣ ≤ C
(
k + h2

)
.

Thus |e1ui| ≤ C(k + h2), i = 1, 2, . . . I − 1.
Likewise, ∣∣e1vi∣∣ ≤ C

(
k + h2

)
, i = 1, 2, . . . , I − 1.



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1529

Now assume that |esui| ≤ Cs

(
k + h2

)
and |esvi| ≤ Cs

(
k + h2

)
, where s = 0, 1, 2, . . . , n,

Cs > 0, C∗ = max0≤s≤nCs.
For n+ 1, we set |en+1

uj | = max1≤i≤I−1 |en+1
ui |, and so∣∣∣en+1

uj

∣∣∣ = (1 + 2rn)
∣∣∣en+1

uj

∣∣∣− rn

(∣∣∣en+1
uj

∣∣∣+ ∣∣∣en+1
uj

∣∣∣)
≤ (1 + 2rn)

∣∣∣en+1
uj

∣∣∣− rn

(∣∣∣en+1
uj+1

∣∣∣+ ∣∣∣en+1
uj−1

∣∣∣)
≤
∣∣∣(1 + 2rn)e

n+1
uj − rn

(
en+1
uj+1 + en+1

uj−1

)∣∣∣
=
∣∣∣enuj + kn

(
F
(
δxu

n
j , v

n
j

)
− F

(
δxU

n
j , V

n
j

))
+ Tn+1

j

∣∣∣ .
Thus ∣∣∣en+1

uj

∣∣∣ ≤ ∣∣enuj∣∣+ kn
∣∣F (δxunj , vnj )− F

(
δxU

n
j , V

n
j

)∣∣+ ∣∣∣Tn+1
j

∣∣∣ .
By Lemma (1), we obtain∣∣∣en+1

uj

∣∣∣ ≤ ∣∣enuj∣∣+ knL1

∣∣δx (unj − Un
j

)∣∣+ knL2

∣∣vnj − V n
j

∣∣+ ∣∣∣Tn+1
j

∣∣∣
≤
∣∣enuj∣∣+ kL1

∣∣δxenuj∣∣+ kL2

∣∣envj∣∣+ ∣∣∣Tn+1
j

∣∣∣ .
By setting L = Max {L1, L2} , we get∥∥En+1

u

∥∥ ≤ ∥En
u∥+ 2kL1 ∥En

u∥+ kL2 ∥En
v ∥+

∣∣∣Tn+1
j

∣∣∣
≤ C∗ (k + h2

)
+ 3kLC∗ (k + h2

)
+ C

(
k + h2

)
≤ (1 + 3kL)C∗ (k + h2

)
+ C

(
k + h2

)
= [(1 + 3kL)C∗ + C]

(
k + h2

)
.

It follows that
∥∥En+1

u

∥∥ ≤ C
(
k + h2

)
, n = 0, 1, . . ., i = 1, 2, . . . , I − 1.

Similarly, we can show that∥∥En+1
v

∥∥ ≤ C
(
k + h2

)
, n = 0, 1, . . . .

Definition 2. If the specified conditions are satisfied, a solution of both the Explicit and
Implicit Euler schemes will together blow up within a finite time denoted as Th:

• ∥Un
h ∥∞ → ∞, ∥V n

h ∥∞ → ∞, as n → ∞,

• Th =
∑∞

n kn.

Remark 1.

• The matrix A = (I − rnhH) exhibits diagonal dominance, characterized by the pres-
ence of positive real diagonal components. This observation suggests that matrix A
is positive-definite and non singular [26], Therefore, it is possible to solve the linear
systems(3.4)-(3.5) and get a unique solution.



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1530

• The finite difference schemes, explicit Euler schemes, and implicit Euler schemes
that have been developed demonstrate consistency, stability, and convergence.

• Explicit (implicit) Euler numerical systems provide approximate solutions for each
set time interval l [0, T ], with a convergence rate of O

(
k + h2

)
where k = maxn kn,

However, when using the time-steps formulas (2.3) and (3.3), the convergence rate
becomes o (hα), as h → 0, for α ≤ 2. The numerical blow-up time is likely to exhibit
a similar pattern of convergence.

• The numerical blow-up time of a fully discrete formula, such as explicit Euler or
implicit Euler, is considered the blow-up time for system (1.1).

• The duration of numerical blow-up achieved through the utilization of a discrete
scheme is contingent upon both the space-step h and the selection of time-steps
kn.Nevertheless, it should be noted that the time-stepping formulas (2.3) and (3.3)
are also contingent upon the space-step. Hence, in the examination of the numerical
convergence of the aforementioned schemes, namely explicit Euler and implicit Eu-
ler, it is sufficient to enhance the space-stepping technique in order to calculate the
error boundaries and numerical order of convergence.

4. Numerical experiments

In this section, we present an estimation of the numerical blow-up times for two numer-
ical experiments, employing distinct space steps, utilizing the finite difference techniques
described in this study, namely explicit Euler and implicit Euler. The Matlab (R2020a)
software is used to write all the numerical computing codes.We use h to denote the space-
step, and kn. is the time-step, h = 1

I . Furthermore, we assess the rate at which the
numerical blow-up time increases for each of these strategies. The numerical blow-up time
is determined by identifying a value of m ∈ N such that the condition:∥Um

h ∥∞ ≥ 1015 or
∥V n

h ∥∞ ≥ 1015 holds . Furthermore, the numerical blow-up time to the examined problem
is denoted as Th = tm =

∑m
n=0 kn The error bounds between T2h and Th are denoted as

. Eh = |T2h − Th|, .The numerical results obtained using the suggested schemes (explicit
Euler and implicit Euler) are presented in tables for each case. The mesh sizes used are
I = {20, 40, 80, 160, 320}, The tables display the iteration count, occurrence of numeri-
cal blow-up, numerical blow-up durations, central processing unit (CPUT) durations in
seconds, and the error limits of numerical blow-up durations. In order to empirically inves-
tigate the rate of numerical convergence for numerical blow-up durations, it is necessary to
consider various mesh sizes with specific values ofα,The formula employed for this purpose
is Sh = log(E2h/Eh)

log 2 . In addition, numerical simulations are conducted to corroborate the
numerical findings.



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1531

4.1. Numerical Experiments

Example 1.

ut = uxx + v5 − |ux|1.6 , vt = vxx + u6 − |vx|1.6 , x ∈ (0, 1), t ∈ (0, T ),
u(0, t) = u(1, t) = 0,
v(0, t) = v(1, t) = 0, t ∈ (0, T ),
u(x, 0) = 80

(
x− x2

)
, v(x, 0) = 90

(
x− x2

)
, x ∈ (0, 1).

 (4.1)

Example 2.

ut = uxx + v3 − |ux|1.1 , vt = vxx + u4 − |vx|1.2 , x ∈ (0, 1), t ∈ (0, T ),
u(0, t) = u(1, t) = 0,
v(0, t) = v(1, t) = 0, t ∈ (0, T ),
u(x, 0) = 100(sinπx), v(x, 0) = 90(sinπx), x ∈ (0, 1).

 (4.2)



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1532

Table 1: Example 1, Explicit formula, α = 1

h m Th CPUT Eh Sh

1/20 3 6.81E-04 0.088252 · · · · · ·
1/40 4 2.09E-04 0.082713 4.72E-04 · · ·
1/80 4 5.29E-05 0.127637 1.56E-04 1.5971
1/160 4 1.42E-05 0.154282 3.87E-05 2.0095
1/320 4 4.73E-06 0.261390 9.44E-06 2.0372

Table 2: Example 1, Explicit formula, α = 2

h m Th CPUT Eh Sh

1/20 4 3.45E-05 0.096988 · · · · · ·
1/40 4 6.49E-06 0.081515 2.80E-05 · · ·
1/80 5 1.44E-06 0.104819 5.05E-06 2.4727
1/160 7 5.16E-07 0.151543 9.26E-07 2.4473
1/320 15 2.88E-07 0.486987 2.28E-07 2.0231

Table 3: Example 1, Implicit formula, α = 1

h m Th CPUT Eh Sh

1/20 3 6.81E-04 0.333453 · · · · · ·
1/40 4 2.09E-04 0.242359 4.72E-04 · · ·
1/80 4 5.29E-05 0.447964 1.56E-04 1.5972
1/160 4 1.42E-05 1.003357 3.87E-05 2.0094
1/320 5 4.73E-06 2.699985 9.44E-06 2.0373



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1533

Table 4: Example 1, Implicit formula, α = 2

h m Th CPUT Eh Sh

1/20 4 3.45E-05 0.202292 · · · · · ·
1/40 4 6.49E-06 0.219009 2.80E-05 · · ·
1/80 5 1.44E-06 0.259753 5.05E-06 2.4727
1/160 7 5.16E-07 0.847733 9.26E-07 2.4473
1/320 15 2.88E-07 2.355910 2.28E-07 2.0231

Table 5: Example 2, Explicit formula, α = 1

h m Th CPUT Eh Sh

1/20 4 1.59E-04 0.337957 · · · · · ·
1/40 5 5.73E-05 0.086509 1.02E-04 · · ·
1/80 6 2.135E-05 0.135266 3.59E-05 1.5083
1/160 7 9.15E-06 0.165087 1.22E-05 1.5590
1/320 9 5.01E-06 0.316164 4.14E-06 1.5566

Table 6: Example 2, Explicit formula, α = 2

h m Th CPUT Eh Sh

1/20 7 1.01E-05 0.785978 · · · · · ·
1/40 12 4.15E-06 0.102545 5.96E-06 · · ·
1/80 36 3.145E-06 0.113383 1.01E-06 2.5626
1/160 163 2.93E-06 0.155124 2.17E-07 2.2135
1/320 879 2.89E-06 0.424556 4.31E-08 2.3346

4.2. Discussion

By the numerical results in Examples 1 and 2, we refer to the following comments:

• The numerical blow-up can happen simultaneously at one time and one point (
x = 0.5), and that agreed with the known theoretical blow-up results of the single
equation of system (1.1), [25].

• When the space-steps are refined, the blow-up time errors-bounds decrease. This
refers that the numerical blow-up times sequence Th, is convergent, as the space-step
approaches to 0



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1534

Table 7: Example 2, Implicit formula, α = 1

h m Th CPUT Eh Sh

1/20 4 1.59E-04 0.360841 · · · · · ·
1/40 5 5.73E-05 0.237088 1.02E-04 · · ·
1/80 6 2.13E-05 0.329224 3.59E-05 1.5084
1/160 7 9.15E-06 0.855710 1.22E-05 1.5589
1/320 9 5.01E-06 2.501314 4.14E-06 1.5567

Table 8: Example 2, Implicit formula, α = 2

h m Th CPUT Eh Sh

1/20 6 1.01E-05 0.459823 · · · · · ·
1/40 12 4.15E-06 0.246959 5.96E-06 · · ·
1/80 36 3.15E-06 0.305770 1.01E-06 2.5627
1/160 163 2.93E-06 0.858916 2.17E-07 2.2136
1/320 879 2.89E-06 2.554193 4.31E-08 2.3346

Figure 1: development of numerical blow-up solution with passage of time emerging from using Explicit for
Example 1, withh = 320, α = 2

• Convergence order of the numerical blow-up times, Sh is close to or larger than the
value of α, this refers, the numerical order of convergence is: 0 (hα+ϵ), where ϵ > 0

• With the time-stepping formulas (2.3) and (3.3), the number of iteration that is
required to achieve blow-up is increasing as the value of α increases.

• When we refine the spatial step, or compare CPUT of implicit method with that of
explicit method, we find the CPUT times are increasing.



M.I.Khalil et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1516-1538 1535

Figure 2: development of numerical blow-up solution with passage of time emerging from using Implicit for
Example 1, with h = 320, α = 2

Figure 3: development of numerical blow-up solution with passage of time emerging from using Explicit for
Example 2, with h = 320, α = 2

Figure 4: development of numerical blow-up solution with passage of time emerging from using Implicit for
Example 2— with h = 320, α = 2

• From tow figures (1) and (2), we find in each of studied problem, the numerical
blow-up growth-rates, obtained by use explicit Euler scheme, is almost the same as
that obtained by use implicit Euler scheme.



REFERENCES 1536

5. Conclusions

This paper deals with numerical approximations of a one-dimensional coupled reaction-
diffusion system with gradient terms. Namely, we propose both Euler explicit and implicit
finite difference methods with a non-fixed time-stepping procedure to estimate the numer-
ical blow-up time of the considered problem. Moreover, some numerical experiments are
given to illustrate the efficiency, accuracy, and numerical order of convergence of the
proposed technique. The obtained numerical results show that the used finite difference
schemes with the proposed non-fixed time-stepping procedure can give accurate results
with high order of numerical convergence. Furthermore, the numerical results show that
the blow-up can only occur at the center point. The two proposed schemes are effective
and can be easily used comprised with other numerical techniques. As future work, we
may use Crank-Nicolson technique to computer the numerical solution of the considered
system.

Acknowledgements

The authors would like to thank Universiti Kebangsaan Malaysia (UKM) for support-
ing their work.

References

[1] Luis M Abia, JC Lopez-Marcos, and Julia Mart́ınez. Blow-up for semidiscretiza-
tions of reaction-diffusion equations. Applied numerical mathematics, 20(1-2):145–
156, 1996.

[2] Omar Abu Arqub. Numerical simulation of time-fractional partial differential equa-
tions arising in fluid flows via reproducing kernel method. International Journal of
Numerical Methods for Heat & Fluid Flow, 30(11):4711–4733, 2020.

[3] EA Az-Zo’bi. A reliable analytic study for higher-dimensional telegraph equation. J.
Math. Comput. Sci, 18(4):423–429, 2018.

[4] Emad Az-Zo’bi, Ahmet Yildirim, and Lanre Akinyemi. Semi-analytic treatment of
mixed hyperbolic–elliptic cauchy problem modeling three-phase flow in porous media.
International Journal of Modern Physics B, 35(29):2150293, 2021.

[5] M Chipot and FB Weissler. Some blowup results for a nonlinear parabolic equation
with a gradient term. SIAM journal on mathematical analysis, 20(4):886–907, 1989.

[6] Miroslav Chlebik, Marek Fila, and Pavol Quittner. Blow-up of positive solutions of a
semilinear parabolic equation with a gradient term. Dynamics of Continuous, Discrete
and Impulsive Systems Series A: Mathematical Analysis, 10:525–537, 02 2003.

[7] Keng Deng. Blow-up rates for parabolic systems. Zeitschrift für angewandte Mathe-
matik und Physik ZAMP, 47:132–143, 1996.



REFERENCES 1537

[8] Marek Fila. Remarks on blow up for a nonlinear parabolic equation with a gradient
term. Proceedings of the American Mathematical Society, 111(3):795–801, 1991.

[9] Avner Friedman and Yoshikazu Giga. A single point blow-up for solutions of semi-
linear parabolic systems. J. Fac. Sci. Univ. Tokyo Sect. IA Math, 34(1):65–79, 1987.

[10] Hiroshi Fujita. On the blowing up of solutions of the Cauchy problem for ut =
δu+ u1+α. J. Fac. Sci. Univ. Tokyo Sect. I, 13:109, 1966.

[11] Houda Hani and Moez Khenissi. On a finite difference scheme for blow up solutions
for the Chipot–Weissler equation. Applied Mathematics and Computation, 268:1199–
1216, 2015.

[12] Houda Hani and Moez Khenissi. Blow-up of semi-discrete solution of a nonlinear
parabolic equation with gradient term. arXiv preprint arXiv:2010.08867, 2020.

[13] Manar Ismael, Ishak Hashim, Maan Rasheed, and Eddie Ismail. Numerical finite
difference approximations of a coupled parabolic system with blow-up. Journal of
Mathematics and Computer Science, 32:387–407, 11 2023.

[14] Stanley Kaplan. On the growth of solutions of quasi-linear parabolic equations. Com-
munications on Pure and Applied Mathematics, 16(3):305–330, 1963.

[15] Monica Marras, S Vernier-Piro, and Giuseppe Viglialoro. Estimates from below of
blow-up time in a parabolic system with gradient term. Int. J. Pure Appl. Math,
93(2):297–306, 2014.

[16] Niraj Mehta, Vipul B Gondaliya, and Jayesh Gundaniya. Applications of differ-
ent numerical methods in heat transfer-a review. International Journal of Emerging
Technology and Advanced Engineering, 3(2):363–368, 2013.

[17] Pavol Quittner and Philippe Souplet. Superlinear parabolic problems. Springer, 2019.

[18] Maan A Rasheed. Blow-up properties of a coupled system of reaction-diffusion equa-
tions. Iraqi journal of Science, pages 3052–3060, 2021.

[19] Maan A Rasheed. On blow-up solutions of a parabolic system coupled in both equa-
tions and boundary conditions. Baghdad Science Journal, 18(2):0315–0315, 2021.

[20] Maan A Rasheed, Hassan Abd Salman Al-Dujaly, Talat Jassim Aldhlki, et al. Blow-
up rate estimates for a system of reaction-diffusion equations with gradient terms.
International Journal of Mathematics and Mathematical Sciences, 2019, 2019.

[21] Maan A Rasheed and Luma J Barghooth. Blow-up set and upper rate estimate for
a semilinear heat equation. In Journal of Physics: Conference Series, volume 1294,
page 032013. IOP Publishing, 2019.



REFERENCES 1538

[22] Maan A Rasheed and Miroslav Chlebik. Blow-up rate estimates and blow-up set for a
system of two heat equations with coupled nonlinear neumann boundary conditions.
Iraqi journal of Science, pages 147–152, 2020.

[23] Maan A Rasheed, Raad A Hameed, Sameer K Obeid, and Ali F Jameel. On numerical
blow-up solutions of semilinear heat equations. Iraqi Journal of Science, 61(8):2077–
2086, 2020.

[24] Maan A Rasheed, Raad Awad Hameed, and Amal Nouman Khalaf. Numerical blow-
up time of a one-dimensional semilinear parabolic equation with a gradient term.
Iraqi Journal of Science, pages 354–364, 2023.

[25] Philippe Souplet. Recent results and open problems on parabolic equations with gra-
dient nonlinearities. Electronic Journal of Differential Equations (EJDE)[electronic
only], 2001:Paper–No, 2001.

[26] Richard S. Varga. Matrix Iterative Analysis. Prentice-Hall Series in Automatic Com-
putation. Prentice-Hall, Englewood Cliffs, 1962.

[27] Pinghui Zhuang and Fawang Liu. Finite difference approximation for two-dimensional
time fractional diffusion equation. Journal of Algorithms & Computational Technol-
ogy, 1(1):1–16, 2007.