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

Finite Difference Scheme for One-Dimensional Coupled
Parabolic System with Blow-up

Manar I. Khalil1,2,∗, Ishak Hashim1, Maan A. Rasheed3, Shaher Momani4

1 Department of Mathematical Sciences, Faculty of Science and Technology,
Universiti Kebangsaan Malaysia, 43600 UKM Bangi, Selangor, Malaysia
2 Department of Applied Earth Sciences, College of Science, Tikrit University, Tikrit, Iraq
3 Department of Mathematics, College of Science for Women, University of Baghdad,
Baghdad, Iraq
4 Nonlinear Dynamics Research Center (NDRC), Ajman University, Ajman,
United Arab Emirates

Abstract. This study aims to find an efficient technique to estimate the blow-up time (BUT)
of one-dimensional semi-linear coupled parabolic systems. Firstly, a fully discrete finite difference
formula is derived with a non-fixed time-stepping formula, based on the Crank-Nicolson method.
In addition, the consistency, stability, and convergence of the proposed scheme are considered.
Secondly, two numerical experiments are presented. For each experiment, we apply the proposed
scheme to calculate the numerical blow-up time, error bounds, and the numerical order of conver-
gence for blow-up times. The obtained results show that the proposed C.N scheme is consistent
with the system considered. However, it is conditionally stable, and the Crank-Nicolson scheme
converges in the stability region and achieves first- and second-order accuracy in temporal and
spatial dimensions, respectively. Also, it helps to increase the order of numerical convergence.
Furthermore, the numerical experiments demonstrate that the numerical blow-up simultaneously
occurs at only the center point. Finally, the numerical blow-up time sequence is convergent. More-
over, convergence for the blow-up time agrees well with the theoretical order of convergence of the
proposed scheme.

2020 Mathematics Subject Classifications: 35K57, 65M06, 65M12, 35B44

Key Words and Phrases: Fully discrete, Crank-Nicolson formula, Numerical blow-up times,
Convergence, Blow-up time, Consistency, Stability

1. Introduction

In recent decades, many real-world problems have been modeled by PDEs, see [1–
3]. There are many parabolic-type PDEs whose solutions cannot possibly be extended

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6140

Email addresses: p119266@siswa.ukm.edu.my (M. I. Khalil)

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 2 of 17

globally over time. This phenomenon is known as blow-up (BU), and it occurs in semi-
linear diffusion equations, when the heat source is strong enough [4–8] this phenomenon
can be applied to interconnected systems, where each variable that relies on another
experiences finite growth within a designated time frame. In this specific situation, it is
known that the dependent variables in the equations demonstrate synchronized increases
[9, 10].This phenomenon has several applications, especially in combustion theory and
heat propagation [11] .
According to [10] , a solution of a diffusion equation defined on Ω × {t > 0} blows up, if
there exist T < ∞, such that it is still bounded as long as 0 < t < T , while it is unbounded
as t is close to the blow-up time T ,

i.e. sup
x∈Ω

|u(x, t)| → ∞ as t → T−.

For a coupled diffusion system:

ut = uxx + F (u, v), vt = vxx +G(u, v), (x, t) ∈ Ω× (0, T ), F,G : R2 → R,

a solution (u, v) blows up simultaneously if there exists T < ∞ such that u and v blow up
at T .
i.e. supx∈Ω |u(x, t)| → ∞, & supx∈Ω |v(x, t)| → ∞, as t → T−,
while for t < T, supx∈Ω{|u(x, t)|, |v(x, t)|} < ∞,

In this work, we consider the system:

ut = uxx + f(v), vt = vxx + g(u), (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)

where u0(x) ∈ C2(R), v0(x) ∈ C2(R), satisfying u0(0) = u0(1) = 0, v0(0) = v0(1) =
0; f, g ∈ C1(R)∩C2(R\{0}), are positive, super-linear differentiable, and increasing func-
tions on (0,∞).

To demonstrate the local existence and distinguish positive solutions to the problem,
some authors used the standard parabolic theory [12]. Moreover, for various nonlinear
functions: f, g, when the initial functions (u0, v0) are sufficiently large, a blow-up may
occur in finite time[9]. In addition, because the system (1) is coupled, only a simultaneous
blow-up occurs. Regarding the blow-up set and blow-up rate estimate, the system (1) has
been extensively studied in [13–19]. Furthermore, imposing some related suggestions in
f, g, it is shown that the (BU) occurs at a single point [15]. Particularly, the following two
instances of f, g are mostly studied [9, 13, 15] :

f(v) = vp, g(u) = uq, p, q > 1
f(v) = epv, g(u) = equ, p, q > 0

Since the last decades, several numerical approximations have been proposed for many
parabolic problems with BU to compute the NBU solution and estimate BUT [20–26].
In [25], the authors proposed an approximate explicit scheme for the system (1) with



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 3 of 17

nonlinear power-type functions, subject to uniform temporal grids. In addition, the com-
putation of the BU time and BU behaviors is studied. Moreover, the convergence of the
NBU time is considered. In addition, the BU set, the BU rate and the BU in the Lp-norm
were taken into account. Furthermore, the relationship between the BU of the numerical
solution and that of the exact solution has been investigated.

Recently, in [26], the authors have studied the semi-discrete problem of system (1)
and its properties regarding convergence and (BUT). Namely, it has been shown that the
(NBUT) and NBUS of a semi-discrete form of (1) converge to the theoretical ones dur-
ing grid refinement. Moreover, they have proposed explicit and implicit finite difference
approximation schemes with non-fixed time-step to estimate the (BUT) to the system
(1).Moreover, they investigated the stability, convergence, and consistency of the sug-
gested methods. In addition, some numerical examples are presented in the form of tables
and figures.

This work aims to find an efficient technique to estimate the (BUT) of the system (1).
Namely, to show the accuracy and efficiency of the C.N. scheme in computing BU solution
and BUT.
This study has six sections. In the second section, some results of the semi-discrete
problem on the problem (1) are recalled. The Crank-Nicolson approximation formulas of
the system (1) are derived in Section three. In the fourth section, consistency, stability,
and convergence are considered. In Section five, two numerical experiments are given
to estimate the (NBUT), error bounds, and numerical order of convergence. Lastly, the
important conclusions and future work are stated in the last section.

1.1. The semi-discrete problem

Let I be a positive integer and xi = ih, 0 ≤ i ≤ I, where h = 1
I , then we can adjective

by the solution: (Uh(t), Vh(t)) ,

Uh(t) = (U0(t), U1(t), . . . , UI(t))
T , Vh(t) = (V0(t), V1(t), . . . , VI(t))

T . (2)

By replacing the second-order space derivative in the problem (1) by the standard second-
order central finite difference operator δ2 [27], we obtain the semi-discrete problem:

d

dt
Ui =

Ui+1 − 2Ui + Ui−1

h2
+ f (Vi) (3)

d

dt
Vi =

Vi+1 − 2Vi + Vi−1

h2
+ g (Ui) (4)

U0(t) = UI(t) = 0, V0(t) = VI(t) = 0, (5)

Ui(0) = U0 (xi) , Vi(0) = V0 (xi) , 0 ≤ i ≤ I. (6)

Definition 1. [26] Let (Uh, Vh) be a nonnegative solution to problem (1)–(2). We say
that (Uh, Vh) achieves blow-up simultaneously in finite time, if there exists Th < ∞ such



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 4 of 17

that:

∥Uh(t)∥∞ < ∞ and ∥Vh(t)∥∞ < ∞ for t ∈ [0, Th) , (7)

∥Uh(t)∥∞ → ∞, ∥Vh(t)∥∞ → ∞ as t → Th, (8)

where ∥Uh(t)∥∞ = max0≤i≤1 |Ui(t)| and ∥7Vh(t)∥∞ = max0≤i≤1 |Vi(t)|.

Theorem 1. [26] Let (Uh, Vh) be a non-negative solution of the semi-discrete problem,
where the reaction functions f, g are locally Lipschitz continuous such that f, g > 0 in
(0,∞). If (Uh, Vh) achieves blow-up at time: Th < ∞, then Th is bounded from below.

Corollary 1. [26] Let (Uh, Vh) be a non-negative solution of the semi-discrete problem
(1)-(2), where f = vp, g = uq, p, q > 1. If (Uh, Vh) achieves blow-up at Th < ∞, then∫ ∞

∥V0∥

ds

G2(s)
=

∫ ∞

∥U0∥

ds

G1(s)
= T ≤ Th, (9)

where

G1(s) =

[
p+ 1

q + 1
Sq+1 − (p+ 1)C0

] p
p+1

, (10)

G2(s) =

[
q + 1

p+ 1
Sp+1 − (q + 1)C0

] q
q+1

, (11)

C0 =
1

q + 1
uq+1
0 − 1

p+ 1
vp+1
0 . (12)

Theorem 2. [26] Assume the following:

a) f, g ∈ C1([0,∞]) are convex (f ′′, g′′ ≥ 0, ) , f(z) > 0 in (0,∞), and ∀ε > 0, we have∫ ∞

ε

dZ

f(Z)
< ∞,

∫ ∞

ε

dZ

g(Z)
< ∞. (13)

b) ∃Z0 ≥ 0 such that(
f(z)
Z

)
≥ λn,

(
g(z)
Z

)
≥ λn for Z ∈ [Z0,∞) , and

limZ→∞
Z

f(Z) = 0, lim Z→∞
Z

g(Z) = 0.

c) The initial conditions are non-negative and such that

Uh(0) ̸= 0, Vh(0) ̸= 0, Q1(0) ≥ Z0, Q2(0) ≥ Z0, (14)

f (Q2(0))− λnQ1(0) > 0, g (Q1(0))− λnQ2(0) ≥ 0, (15)

Where

Q1(t) =
I−1∑
i=1

hUi(t)Υi, Q2(t) =
I−1∑
i=1

hVi(t),Υi, (16)



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 5 of 17

Υi =
sin(πih)∑I−1

m=1 sin(πmh)
,−δ2Υi = λnΥi, i = 1, 2, . . . , I − 1, Υ0 = ΥI = 0,

λn =

(
4

h2

)
sin2

(
πh

2

)
.

(17)

Then the non-negative solution (Uh, Vh) of the discrete problem achieves blow-up at time
Th with

Th ≤
∫ ∞

Q1(0)

dZ1

f (Z2)− λnZ1
, ∀Z2 ∈ [Z0,∞) , (18)

or

Th ≤
∫ ∞

Q2(0)

dZ2

g (Z1)− λnZ2
, ∀Z1 ∈ [Z0,∞) . (19)

Theorem 3. [26] Assume a- The reaction functions: f, g ∈ C1([0,∞]), and problem of
system (1) has a solution: (u, v), u, v ∈ C4,1([0, 1]× [0, T ]).
b- The initial condition

(
U0
h , V

0
h

)
satisfies

∥∥U0
h − uh(0)

∥∥
∞ = O(1), h → 0∥∥V 0

h − vh(0)
∥∥
∞ = O(1), h → 0

(20)

Then, for h sufficiently small, the semi-discrete problem (1) has a unique solution:
(Uh, Vh) , Uh, Vh ∈ C1([0, T ]), RJ+1 such that

maxt∈[o,T ] ∥Uh(t)− uh(t)∥∞ = O
(
∥e(0)∥∞ + h2

)
, h → 0

maxt∈[o,T ] ∥Vh(t)− vh(t)∥∞ = O
(
∥e(0)∥∞ + h2

)
, h → 0

(21)

where ∥e(0)∥∞ = max
{∥∥U0

h − uh(0)
∥∥
∞ ,
∥∥V 0

h − vh(0)
∥∥
∞
}

Theorem 4. [26] Assume
a) The functions f, g ∈ C1([0,∞), R) and f(Z), g(Z) > 0 in (0,∞).
b) There exists Z̄ ≥ 0, such that

•
(
f(Z)
Z

)
≥ π2,

(
g(Z)
Z

)
≥ π2, for Z ∈ [Z̄,∞)

• limZ→∞
Z

f(Z) = 0, limZ→∞
Z

g(Z) = 0

• f(Z̄)− π2Z̄ > 0 , g(Z̄)− π2Z̄ > 0

c) There exists T < ∞ such that, u, v ∈ C4,2([0, 1]× [0, T )) and

lim
t→T

∫ 1

0
u(x, t)w(x)dx = ∞

lim
t→T

∫ 1

0
v(x, t)w(x)dx = ∞

(22)



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 6 of 17

If ∥Uh(0)− uh(0)∥∞ = 0(1), h → 0, then the solution of the discrete problem achieves
blow-up, for h sufficiently small at Th and

lim
h→0

Th = T

2. Crank-Nicolson Method

Define the grid points xi = ih, 0 ≤ i ≤ I, where I ∈ Z+, h = 1/I, tn+1 = tn + kn, kn >
0,∀n.

Set Un
h =

(
Un
1 , U

n
2 . . . Un

I−1

)T
, V n

h =
(
V n
1 , Un

2 . . . V n
I−1

)T
as the numerical solution of problem (1).
In order to derive the Crank-Nicolson method formally, we approximate the partial deriva-

tives in the semi-discrete problem (3),(4), at the mesh points:
(
xi, tn+ 1

2

)
.

Firstly, we approximate the time-derivatives using the 1s order central finite difference
formula:

∂u

∂t

∣∣∣∣n+ 1
2

i

=
1

kn

(
un+1
i − uni

)
+O (Kn) ,

∂v

∂t

∣∣∣∣n+ 1
2

i

=
1

kn

(
vn+1
i − vni

)
+O (kn) ,

and δ2xU(t), δ2xV (t) are approximated as follows:

δ2xUi

(
tn+1/2

)
=

1

2

[
δ2xUi (tn+1) + δ2xUi (tn)

]
+O

(
k2n
)

δ2xVi

(
tn+1/2

)
=

1

2

[
δ2xVi (tn+1) + δ2xVi (tn)

]
+O

(
k2n
)

It follows that

δ2xUi

(
tn+1/2

)
=

1

2h2
[(
Un
i+1 − 2Un

i + Un
i−1

)
+
(
Un+1
i+1 − 2Un+1

i + Un+1
i−1

)]
δ2xVi

(
tn+1/2

)
=

1

2h2
[(
V n
i+1 − 2V n

i + V n
i−1

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

i + vn+1
i−1

)]
While the nonlinear terms are taken as follows:

f
(
V

n+1/2
i

)
= f (V n

i ) +O (Kn) ,

g
(
U

n+1/2
i

)
= g (Un

i ) +O (kn)

By substituting all these formulas in (3),(4), we obtain

1

kn

(
Un+1
i − Un

i

)
=

1

2h2
[(
Un
i+1 − 2Un

i + Un
i−1

)
+
(
Un+1
i+1 − 2Un+1

i + Un+1
i−1

)]
+ f (V n

i )

1

kn

(
V n+1
i − V n

i

)
=

1

2h2
[(
V n
i+1 − 2V n

i + V n
i−1

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

i + V n+1
i−1

)]
+ g (Un

i ) .



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 7 of 17

The last difference equation can be written as follows:

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

2

(
Un+1
i+1 + Un+1

i−1

)
= (1− rn)U

n
i +

rn
2

(
Un
i+1 + Un

i−1

)
+ knf (V n

i ) (23)

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

2

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

i−1

)
= (1− rn)V

n
i +

rn
2

(
V n
i+1 + V n

i−1

)
+ kng (U

n
i ) , (24)

when rn = kn
h2 .

Moreover, one can choose the time step:

kn = min

(
h2,

hα∥∥Un
h

∥∥
∞
,

hα∥∥V n
h

∥∥
∞

)
, 1 ≤ α. (25)

System (1) can be presented in matrix form as follows:(
I +

rn
2
H
)
Un+1
h =

(
I − rn

2
H
)
Un
h + knF

n (26)

(
I +

rn
2
H
)
V n+1
h =

(
I − rn

2
H
)
V n
h + knG

n (27)

where H =

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

. . .
. . .

0 · · · 0 −1 2

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

Fn =
(
f (V n

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

(
V n
I−1

))T
, Gn =

(
g (Un

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

(
Un
I−1,

))T
Lemma 1. Let f, g be differentiable real functions. Then, there exist positive constants
L1, L2, > 0, such that ∣∣∣f (V n

i )− f
(
Ṽ n
i

)∣∣∣ ≤ L1

∣∣∣V n
i − Ṽ n

i

∣∣∣ (28)∣∣∣g (Un
i )− g

(
Ũn
i

)∣∣∣ ≤ L2

∣∣∣Un
i − Ũn

i

∣∣∣ (29)

where Un
i , V

n
i , Ũn

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

Proof. To end this, we use the mean value theorem:∣∣∣f (V n
i )− f

(
Ṽ n
i

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

i − Ṽ n
i

∣∣∣ (30)

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

i .
Sine V n

h , Ṽ n
h are bounded, then there exists C1 > 0 such that |f ′ (Zi)| < C1.

Thus, (28) is valid. Similarly, one can show that (29) holds true.



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 8 of 17

2.1. The algorithm steps for Crank Nicolson method

(i) Input h, U0
h , V

0
h ,p, q, α

(ii) Put n = 0;

(iii) Choose kn according to (25).

(iv) Calculate the numerical vectors: Un+1
h , V n+1

h , using the Crank Nicolson formula (26)
and (27).

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

∥V n
h ∥∞ ≥ 1015

(vi) The numerical blow-up time is tm =
∑m

n=0 kn.

3. Consistency, Stability and Convergence of C. N. scheme

In this section, we aim to analyze the Crank-Nicolson scheme in terms of these three
core properties. We will first examine its consistency by evaluating the local truncation
error. Next, we will explore its stability and conclude the convergence of the method,
ensuring that the numerical solution approximates the exact solution.

Theorem 5. Let (Tn
ui, T

n
vi) be the local truncation error of the Crank-Nicolson formulas

at the grid point
(
xi, tn+ 1

2

)
. Then∣∣∣∣Tn+ 1

2
ui

∣∣∣∣ ≤ C1k + C2h
2,

∣∣∣∣Tn+ 1
2

vi

∣∣∣∣ ≤ C3k + C4h
2 , C1, C2, C3, C4 > 0.

Proof. Replacing the precise solution uni = u
(
xi, tn+ 1

2

)
, vni = v

(
xi, tn+ 1

2

)
in the

Crank- Nicolson, yields that:

T
n+ 1

2
ui =

(
un+1
i − uni

)
− kn

2h2
(
un+1
i+1 − 2un+1

i + un+1
i−1

)
+
(
uni+1 − 2uni + uni−1

)
−knf (vni )

It follows that

T
n+ 1

2
ui = kn

∂un+ 1
2

i

∂t
+O

(
k2n
)− kn

[
∂2un+1

i

∂x2
+O

(
k2n + h2

)]

−kn

[
f

(
v
n+ 1

2
i

)
+O (kn)

]
Thus, there exists C > 0 , such that∣∣∣∣Tn+ 1

2
ui

∣∣∣∣ ≤ C
(
k + h2

)
, where k = max

n∈N
kn,



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 9 of 17

By the same way, one can get

∣∣∣∣Tn+ 1
2

vi

∣∣∣∣ ≤ C
(
k + h2

)
Clearly, T

n+ 1
2

ui and T
n+ 1

2
vi approach zero, as the space (time) steps go to zero. Hence, the

C.N. formulas are consistent.

Theorem 6. The Crank-Nicolson scheme is stable if (1− r) ≥ 0, where k = maxn∈N kn,

r =
k

h2

Proof. Let En
u = (enui, i = 1, 2, . . . I − 1) , En

v = (envi, i = 1, 2, . . . I − 1)
where enui = uni −Un

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

solution of problem (1)

Now, For n = 1,
∥∥E1

u

∥∥ = max1≤i≤I

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

∥∥ = max1≤i≤I

∣∣e1vi∣∣.
Clearly, ∣∣e1uj∣∣ = (1 + r0)

∣∣e1uj∣∣− r0
2

(∣∣e1uj∣∣+ ∣∣e1uj∣∣)
≤ (1 + r0)

∣∣e1uj∣∣− r0
2

(∣∣e1uj+1

∣∣+ ∣∣e1uj−1

∣∣)
≤
∣∣∣(1 + r0) e

1
uj −

r0
2

(
e1uj+1 + e1uj−1

)∣∣∣
=
∣∣∣(1− r0) e

0
uj +

r0
2

(
e0uj+1 + e0uj−1

)
+ k0

(
f
(
v0j
)
− f

(
V 0
j

))∣∣∣ .
Since (1− rn) ≥ 0, it follows that∥∥E1

u

∥∥ ≤ (1− r0)
∥∥E0

u

∥∥+ r0
∥∥E0

u

∥∥+ L1k0
∣∣v0j − V 0

j

∣∣∥∥E1
u

∥∥ ≤
∥∥E0

u

∥∥+ Lk
∥∥E0

v

∥∥
≤ (1 + kL)max

{∥∥E0
u

∥∥ , ∥∥E0
v

∥∥} .
By the same way, one can get:∥∥E1

v

∥∥ ≤
∥∥E0

v

∥∥+ kL
∥∥E0

u

∥∥
≤ (1 + kL)max

{∥∥E0
u

∥∥ , ∥∥E0
v

∥∥}
where , L = max {L1, L2}.

Assume that:

∥Es
u∥ ≤ (1 + kL)smax

{∥∥E0
u

∥∥ , ∥∥E0
v

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

v∥ ≤ (1 + kL)smax
{∥∥E0

u

∥∥ ,∥∥E0
v

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

Regarding n+ 1,

set
∥∥En+1

u

∥∥ = max1≤i≤I

∣∣∣en+1
uj

∣∣∣ = ∣∣∣en+1
uj

∣∣∣ , ∥∥En+1
v

∥∥ = max1≤i≤I

∣∣∣en+1
vj

∣∣∣



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 10 of 17

It is clear that ∣∣∣en+1
uj

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

uj

∣∣∣− rn
2

(∣∣∣en+1
uj

∣∣∣+ ∣∣∣en+1
uj

∣∣∣)
≤ (1 + r0)

∣∣∣en+1
uj

∣∣∣− r0
2

(∣∣∣en+1
uj+1

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

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

n+1
uj − rn

2

(
en+1
uj+1 + en+1

uj−1

)∣∣∣
=
∣∣∣(1− rn) e

n
uj +

rn
2

(
enuj+1 + enuj−1

)
+ kn

(
f
(
vnj
)
− f

(
V n
j

))∣∣∣ .
Since (1− rn) ≥ 0, it follows that∥∥En+1

u

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

u∥+ L1kn
∣∣vnj − V n

j

∣∣∥∥En+1
u

∥∥ ≤ ∥En
u∥+ Lk ∥En

v ∥
≤ (1 + kL)n+1max

{∥∥E0
u

∥∥ ,∥∥E0
v

∥∥}
≤ exp((n+ 1)kL)max

{∥∥E0
u

∥∥ ,∥∥E0
v

∥∥}
Thus ∥∥En+1

u

∥∥ ≤ exp (tn+1L)max
{∥∥E0

u

∥∥ ,∥∥E0
v

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

v

∥∥ ≤ exp (tn+1L)max
{∥∥E0

u

∥∥ ,∥∥E0
v

∥∥} .
Based on the stability definition [22] the Crank-Nicolson scheme is stable provided that
r ≤ 1

Theorem 7. Under the condition of theorem (2), the Crank-Nicolson scheme converges
with the order O

(
k + h2

)
, where k = maxn∈N kn.

Proof. Let En
u = (enui, i = 1, 2, . . . I − 1) , En

v = (envi, i = 1, 2, . . . I − 1)
where enui = uni − Un

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

solution of problem (1)
Suppose that e0ui = 0, e0vi = 0 ∀i = 0, 1, . . ..

To prove this theorem, the following inequalities should be held

en+1
uj ≤ C

(
k + h2

)
, en+1

vj ≤ C
(
k + h2

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

Now, For n = 1, set
∣∣∣e1uj∣∣∣ = max1≤i≤I

∣∣e1ui∣∣.
It clear that ∣∣e1uj∣∣ = (1 + r0)

∣∣e1uj∣∣− r0
2

(∣∣e1uj∣∣+ ∣∣e1uj∣∣)
≤ (1 + r0)

∣∣e1uj∣∣− r0
2

(∣∣e1uj+1

∣∣+ ∣∣e1uj−1

∣∣)
≤
∣∣∣(1 + r0) e

1
uj −

r0
2

(
e1uj+1 + e1uj−1

)∣∣∣
=
∣∣∣(1− r0) e

0
uj +

r0
2

(
e0uj+1 + e0uj−1

)
+ k0

(
f
(
v0j
)
− f

(
V 0
j

))
+ T 0

ui

∣∣∣ .



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 11 of 17

It follows that
∣∣∣e1uj∣∣∣ ≤ Lk

∣∣∣e0uj∣∣∣+ ∣∣T 0
ui

∣∣ = ∣∣T 0
ui

∣∣ ≤ C
(
k + h2

)
,

where L = max {L1, L2}.
Hence

∣∣∣e1uj∣∣∣ ≤ C1

(
k + h2

)
, i = 0, 1, . . . , I − 1, C1 = C.

Suppose that ∥Es
u∥ ≤ Cs

(
k + h2

)
∥Es

v∥ ≤ Cs

(
k + h2

)
s = 0, 1, 2 . . . , n Cs > 0.

For n+ 1, set ∥En
u∥ =

∣∣∣en+1
uj

∣∣∣ = max1≤i≤I

∣∣en+1
ui

∣∣
∣∣∣en+1

uj

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

uj

∣∣∣− rn
2

(∣∣∣en+1
uj

∣∣∣+ ∣∣∣en+1
uj

∣∣∣)
≤ (1 + rn)

∣∣∣en+1
uj

∣∣∣− rn
2

(∣∣∣en+1
uj+1

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

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

n+1
uj − rn

2

(
en+1
uj+1 + en+1

uj−1

)∣∣∣
=
∣∣∣(1− rn) e

n
uj +

rn
2

(
enuj+1 + enuj−1

)
+ kn

(
f
(
vnj
)
− f

(
V n
j

))
+ T

n+1/2
ui

∣∣∣ .
Since (1− rn) ≥ 0, it follows that∥∥En+1

u

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

u∥+ L1kn
∣∣vnj − V n

j

∣∣+ ∣∣∣Tn+1/2
ui

∣∣∣∥∥En+1
u

∥∥ ≤ ∥En
u∥+ Lk ∥En

v ∥+
∣∣∣Tn+1/2

ui

∣∣∣
≤ Cn

(
k + h2

)
+ kLCn

(
k + h2

)
+ C

(
k + h2

)
≤ (1 + kL)Cn

(
k + h2

)
+ C

(
k + h2

)
= [(1 + kL)Cn + C]

(
k + h2

)
.

It follows that
∥∥En+1

u

∥∥ ≤ Cn+1

(
k + h2

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

In the same way, one can get∥∥En+1
v

∥∥ ≤ Cn+1

(
k + h2

)
, n = 0, 1, . . .

Remark

• Clearly, the matrix
(
I + rn

2 H
)
is diagonally dominated; hence, it is non-singular. It

follows that the linear systems (26), (27) are uniquely solvable [28]

• The value of (NBUT) is dependent on the space and time steps.

• As it was proved in theorem (3), the Crank-Nicolson numerical scheme gives approx-
imate solutions of the order of convergence O

(
k + h2

)
; k = maxn kn, while, with the

formula of time steps (25), the convergence rate becomes O (hα), as the space step
approaches zero for 1 ≤ α ≤ 2. A similar convergence order is expected for (NBUT).



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 12 of 17

4. Numerical Experiments

In this section, the Crank-Nicolson scheme is used for two numerical experiments,
with various grid sizes: I = {20, 40, 80, 160, 320}, taking α = 1, 2. In addition, Matlab
(R2020a) software is used for writing all the numerical computations codes. For each
example, some tables and figures are presented to show the numerical results obtained by
using the proposed scheme. Each table includes:

• For the first m ∈ N such that the condition ∥Um
h ∥∞ ≥ 1015 and ∥V m

h ∥∞ ≥ 1015

holds, the value Th = tm =
∑m

n=0 kn is considered the (NBUT) for the studied
problems.

• Eh = |T2h − Th| is the error bond between T2h and Th.

• The order of numerical convergence to the (BUT), using the formula: Sh = log (E2h/Eh) / log 2.

• The unit time of the central processor (CPUT) in seconds.

4.1. Examples

Example 1.

ut = uxx + v5, vt = vxx + u6, 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) = 70

(
x− x2

)
, v(x, 0) = 80

(
x− x2

)
, x ∈ (0, 1)

 (31)

Example 2.

ut = uxx + v6 vt = vxx + u7, 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) = 60
(
1− ex

2−x
)
, v(x, 0) = 70

(
1− ex

2−x
)
, x ∈ (0, 1)

 (32)



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 13 of 17

Table 1: Example 1, C.N. scheme: α = 1

h m Th CPUT Eh Sh

1/20 3 7.6650E-04 0.362506 · · · · · ·
1/40 4 2.0966E-04 0.288462 5.5684E-04 · · ·
1/80 4 5.3939E-05 0.346336 1.5572E-04 1.8383
1/160 4 1.5557E-05 0.877538 3.8382E-05 2.0205
1/320 4 6.4259E-06 2.608468 9.1311E-06 2.0716

Table 2: Example 1, C.N. scheme: α = 2

h m Th CPUT Eh Sh

1/20 4 3.9272E-05 0.230858 · · · · · ·
1/40 4 7.6677E-06 0.335248 3.1604E-05 · · ·
1/80 5 1.9252E-06 0.326475 5.7425E-06 2.4604
1/160 8 7.8273E-07 0.789320 1.1425E-06 2.3295
1/320 21 5.2762E-07 2.642523 2.5511E-07 2.1630

Table 3: Example 2, C.N. scheme: α = 1

h m Th CPUT Eh Sh

1/20 3 8.3371E-04 0.401556 · · · · · ·
1/40 3 2.0886E-04 0.395284 6.2485E-04 · · ·
1/80 3 5.2830E-05 0.327884 1.5603E-04 2.0017
1/160 4 1.4065E-05 0.953267 3.8765E-05 2.0090
1/320 4 4.6714E-06 2.653295 9.3936E-06 2.0450



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 14 of 17

Table 4: Example 2, C.N. scheme: α = 2

h m Th CPUT Eh Sh

1/20 4 4.8891E-05 0.223913 · · · · · ·
1/40 4 8.8974E-06 0.233469 3.9994E-05 · · ·
1/80 4 1.7888E-06 0.329078 7.1086E-06 2.4921
1/160 5 4.4920E-07 1.044031 1.3396E-06 2.4078
1/320 7 1.7424E-07 2.458527 2.7496E-07 2.2845

Figure 1: For example 1, (NBU) solution evolution over time using C.N. scheme, with h = 320, α = 2.

Figure 2: For example 2, (NBU) solution evolution over time using C.N. scheme, with h = 320, α = 2.

4.2. Results and Discussion

In Tables 1-4 and Figures 1-2, the following observations can be made: Firstly,the
(NBU) simultaneously occurs near the center point (x = 0.5), and this is consistent with
the known (BU) results of the problem (1), see [15]. In addition, when the space steps
are refined, the (BUT) error bounds decrease. In other words, the (NBUT) sequence Th is
convergent, since the space step is close to zero. Furthermore, the numerical convergence
order of (NBUT) is O (hα+ϵ), where ϵ > 0.Also,a large value of α requires many iterations
to achieve (BU) compared to a small value. Finally, when spatial steps are refined, one
can see an increase in CPUT times.



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 15 of 17

5. Conclusions

This study employs the Crank-Nicolson scheme with a non-fixed time-stepping for-
mula to identify (NBU) solutions and estimate the (NBUT) for a system of two coupled
semi-linear heat equations subject to homogeneous Dirichlet boundary conditions. The
consistency, stability, and convergence properties of the Crank-Nicolson scheme are rigor-
ously examined. To validate the theoretical findings and determine the convergence rate of
the NBUT, two numerical experiments are conducted. The numerical results, presented
through tables and figures, demonstrate that the proposed method achieves a high de-
gree of convergence and exhibits good computational efficiency. Future research directions
may include exploring special approximations for nonlinear terms in the system to enhance
the accuracy of the scheme, potentially achieving second-order truncation errors in both
temporal and spatial dimensions.

References

[1] Maan A Rasheed, Sean Laverty, and Brittany Bannish. Numerical solutions of a linear
age-structured population model. In AIP Conference Proceedings, volume 2096. AIP
Publishing, 2019.

[2] Raheam Al-Saphory, Zinah Khalid, and Abdelhaq EL Jai. Regional boundary gradi-
ent closed loop control system and γ* agfo-observer. In Journal of Physics: Confer-
ence Series, volume 1664, page 012061. IOP Publishing, 2020.

[3] Soraya Rekkab, Samir Benhadid, and Raheam Al-Saphory. An asymptotic analysis of
the gradient remediability problem for disturbed distributed linear systems. Baghdad
Science Journal, 19(6 (Suppl.)):1623–1623, 2022.

[4] Avner Friedman and Bryce McLeod. Blow-up of positive solutions of semilinear heat
equations. Indiana University Mathematics Journal, 34(2):425–447, 1985.

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

[6] Yuzhu Han. Blow-up at infinity of solutions to a semilinear heat equation with loga-
rithmic nonlinearity. Journal of Mathematical Analysis and Applications, 474(1):513–
517, 2019.

[7] Maan A Rasheed, Raad Awad Hameed, Amal Nouman Khalaf, and Hiba Ali Ka-
reem. Estimating the blow-up time for chipot-weissler equation using crank-nicolson
method. J. Appl. Math. & Informatics Vol, 43(1):123–137, 2025.

[8] Manar Khalil, Ishak Hashim, Maan Rasheed, Faieza Samat, and Shaher Momani.
Numerical finite-difference approximations of a coupled reaction-diffusion system with
gradient terms. European Journal of Pure and Applied Mathematics, 17(3):1516–1538,
2024.

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

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



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 16 of 17

[11] A.A. Lacey. Diffusion models with blow-up. Journal of Computational and Applied
Mathematics, 97(1):39–49, 1998.

[12] Olga A Ladyženskaja, Vsevolod Alekseevich Solonnikov, and Nina N Ural’ceva. Linear
and quasi-linear equations of parabolic type, volume 23. American Mathematical Soc.,
1988.

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

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

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

[16] V Galaktionov. Parabolic system of quasilinear equations. i. DIFFERENTIAL
EQUAT., 19(12):1558–1574, 1984.

[17] ZG Lin, CH Xie, and MX Wang. The blow-up rate of positive solutions of a parabolic
system. Northeast. Math. J, 13:327–378, 1997.

[18] Philippe Souplet. Single-point blow-up for a semilinear parabolic system. Journal of
the European Mathematical Society, 11(1):169–188, 2009.

[19] Sining Zheng, Lizhong Zhao, and Feng Chen. Blow-up rates in a parabolic system of
ignition model. Nonlinear Analysis: Theory, Methods & Applications, 51(4):663–672,
2002.

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

[21] Maan A Rasheed and Faez N Ghaffoori. Numerical blow-up time and growth rate
of a reaction-diffusion equation. Italian journal of pure and applied mathematics,
44:805–813, 2020.

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

[23] Weizhang Huang, Jingtang Ma, and Robert D Russell. A study of moving mesh pde
methods for numerical simulation of blowup in reaction diffusion equations. Journal
of Computational Physics, 227(13):6532–6552, 2008.

[24] Chien-Hong Cho and Hisashi Okamoto. Finite difference schemes for an axisym-
metric nonlinear heat equation with blow-up. Electronic Transactions on Numerical
Analysis, 52(2):391–415, 2020.

[25] Chien-Hong Cho and Ying-Jung Lu. On the numerical solutions for a parabolic system
with blow-up. AIMS Mathematics, 6(11):11749–11777, 2021.

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

[27] Gordon D Smith. Numerical solution of partial differential equations: finite difference
methods. Oxford university press, 1985.



M. I. Khalil et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6140 17 of 17

[28] Richard Varga. S. 1962 matrix iterative analysis. Englewood Cliffs: Prenticeflall,
1962.