An Efficient Explicit Scheme for Solving the 2D Heat Equation with Stability and Convergence Analysis

Md Nazmul Haque¹, Ruma Akter², Md Shahadat Hossain Mojumder^3*

¹Department of Mathematical Sciences, Kent State University, Kent, OH, USA.
²Department of Mathematics, Comilla University, Cumilla, Bangladesh.
³School of Computer Science and Technology, University of Bedfordshire, Luton, UK.

DOI: 10.4236/jamp.2025.137127                

Abstract

This paper presents a comprehensive numerical study of the two-dimensional time-dependent heat conduction equation using the Forward Time Centered Space (FTCS) finite difference scheme. The heat equation is a fundamental parabolic partial differential equation, models the diffusion of thermal energy in a medium and is applicable in areas such as thermal insulation design, microchip cooling, and biological heat transfer. Due to the limitations of analytical methods in handling complex geometries and boundary conditions, we employ the FTCS scheme. The problem is formulated with Dirichlet boundary conditions and a sinusoidal initial condition for which an exact analytical solution is known. We derive the FTCS discretization using Taylor series-based approximations and perform a detailed von Neumann stability analysis to establish the Courant-Friedrichs-Lewy (CFL) condition. The scheme’s performance is evaluated through numerical simulations on a uniform grid, with results compared against the exact solution. Simulation results show that the FTCS scheme achieves L² and max-norm errors on the order of 10⁻¹¹ and 10⁻¹⁰, respectively, under stable conditions. Graphical comparisons further demonstrate excellent agreement between numerical and analytical solutions. Overall, the FTCS method proves to be a robust and reliable tool for solving heat conduction problems, provided the stability criterion is satisfied.

Keywords

Explicit Scheme, Finite Difference Method, CFL Condition, Stability, Convergence, Von-Neumann Method, Error Norms

1. Introduction

Heat is a form of energy that transfers from one medium to another due to a temperature difference. According to the second law of thermodynamics, heat naturally flows from a hotter region to a cooler one, proportional to the temperature gradient and the material’s thermal conductivity [1]. This phenomenon governs many real-world applications, from industrial processes to biological systems.

The mathematical modeling of this process is described by the heat equation, a parabolic, second-order linear partial differential equation (PDE) that captures the temporal and spatial evolution of temperature in a medium. It explains how thermal energy diffuses through a region over time. For example, when a hot object is submerged in cooler water, its temperature gradually decreases, and eventually, thermal equilibrium is reached in the absence of external heat sources.

While analytical methods such as separation of variables [2]-[6], Fourier series [7] [8], and domain decomposition [9]-[11] have been effectively used to solve the heat equation under simplified conditions, they are often impractical for complex geometries or boundary conditions. To overcome these limitations, numerical methods, particularly the Finite Difference Method (FDM) [12]-[15], have become widely adopted. FDM approximates the derivatives in the PDE using finite difference expressions derived from Taylor series expansions [16]-[18], converting the PDE into a system of linear algebraic equations that can be efficiently solved using iterative algorithms.

The Forward Time Centered Space (FTCS) scheme is a classical explicit finite difference method known for its simplicity and ease of implementation. Despite its conditional stability, the FTCS scheme remains popular for short-time simulations where the time step can be carefully controlled.

In this paper, we investigate the two-dimensional heat equation using the FTCS method. We derive the numerical scheme, perform a rigorous von Neumann stability analysis, and validate the method by comparing the numerical solution with the exact analytical solution. The error is evaluated both graphically and numerically using maximum and L²-norms. The numerical experiments are carried out in MATLAB, and the results demonstrate the accuracy and reliability of the method under stable conditions.

The remainder of the paper is organized as follows: Section 2 presents the mathematical model and initial-boundary value problem. Section 3 describes the finite difference discretization and development of the FTCS scheme. The stability analysis using von Neumann’s method is given in Section 4. Section 5 provides a detailed error analysis and discussion of the numerical results. Finally, Section 6 offers concluding remarks.

2. Model Equations

In our model problem, we consider the two-dimensional time dependent heat conduction equation of the form [18],

$\frac{\partial u}{\partial t}={\alpha }^2\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right),0<x<a,0<y<b,t>0,$ (1)

with initial condition

$u\left(x,y,0\right)=f\left(x,y\right),0<x<a,0<y<b,$ (2)

and boundary conditions

$u\left(0,y,t\right)=u\left(a,y,t\right)=0,0<y<b,t>0,$ (3)

$u\left(x,0,t\right)=u\left(x,b,t\right)=0,0<x<a,t>0,$ (4)

where, $u=u\left(x,y,t\right)$ is the temperature of the substance, ${\alpha }^2$ the thermal diffusivity of the substance. Equations (1), (2), (3) and (4) together are called initial boundary-value problems.

3. Finite Difference Discretization

The finite difference method (FDM) is a numerical technique that approximates derivatives by using finite differences at discrete grid points. Through discretization, the continuous domain in space and time is converted into a finite grid, allowing partial diffessrential equations like the heat equation to be solved as systems of algebraic equations [4] [14] [15].

In order to describe the finite difference techniques, let us consider the rectangle $\text{Ω}\times \text{Γ}$ into a finite number of nodes $\left(x_i,t_j\right)$ , where Ω denotes the spatial domain and Γ represents the time domain.

3.1. Approximation of the Derivatives of a Function by Finite Differences

In this section, we derive the finite difference approximations for the partial derivatives in the heat Equations (1)-(4).

Step 1: (Discretization of the domain) We divide the spatial domain $\left(x,y\right)$ into a uniform grid with step sizes $\text{Δ}x$ and $\text{Δ}y$ , and we divide the time domain into discrete steps of size $\text{Δ}t$ .

Let $u_{i,j}^n$ denote the temperature at position $\left(x_i,y_j\right)$ at time $t_n$ . The spatial points are indexed as:

$x_i=i\text{Δ}x,i=0,1,2,\cdots ,N_x$

$y_j=j\text{Δ}y,j=0,1,2,\cdots ,N_y$

$t_n=n\text{Δ}t,n=0,1,2,\cdots$

where, $N_x$ and $N_y$ are the number of grid intervals in $x$ and $y$ directions, respectively, and the time step size is defined as $\text{Δ}t=\frac{T}{N_t}$ , where $N_t$ is the total number of time steps.

Step 2: (Finite difference approximation): We use a forward difference scheme for the time derivative (explicit in time) and central difference schemes for the spatial derivatives. These approximations are given as follows:

The first-order time derivative is approximated by the forward difference:

$\frac{\partial u}{\partial t}\approx \frac{u_{i,j}^{n+1}-u_{i,j}^n}{\text{Δ}t}.$ (5)

The second-order spatial derivatives are approximated by the central difference scheme:

$\frac{{\partial }^{2}u}{\partial x^2}\approx \frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\text{Δ}{x}^2},$ (6)

$\frac{{\partial }^{2}u}{\partial y^2}\approx \frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\text{Δ}{y}^2}.$ (7)

3.2. Forward Time Centered Space (FTCS) Explicit Scheme

To construct the FTCS formula for our Equation (1) [6] [14] [15], we substitute the approximations (5)-(7) in (1). We obtain

$\frac{u_{i,j}^{n+1}-u_{i,j}^n}{\text{Δ}t}={\alpha }^2\left(\frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\text{Δ}{x}^2}+\frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\text{Δ}{y}^2}\right).$

Rearranging for $u_{i,j}^{n+1}$ :

$u_{i,j}^{n+1}=u_{i,j}^n+{\alpha }^2\text{Δ}t\left(\frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\text{Δ}{x}^2}+\frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\text{Δ}{y}^2}\right).$

Define stability parameters:

$S_x=\frac{{\alpha }^2\text{Δ}t}{\text{Δ}{x}^2},S_y=\frac{{\alpha }^2\text{Δ}t}{\text{Δ}{y}^2}.$

The update formula simplifies to:

$u_{i,j}^{n+1}=u_{i,j}^n+S_x\left(u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n\right)+S_y\left(u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n\right).$ (8)

Figure 1. FTCS Grid: The node $\left(i,j\right)$ and its four immediate neighbors used in the Forward Time Centered Space (FTCS) discretization for the 2D heat equation.

Figure 1 illustrates the finite difference stencil used in the Forward Time Centered Space (FTCS) scheme for the two-dimensional heat equation. The central node $\left(i,j\right)$ corresponds to the grid point where the temperature is being updated. Its four immediate neighbors $\left(i,j-1\right)$ , $\left(i,j+1\right)$ , $\left(i-1,j\right)$ , and $\left(i+1,j\right)$ are used to approximate the second-order spatial derivatives in the $x$ and $y$ directions via the central difference method. This five-point stencil structure is fundamental to the explicit FTCS update formula in Equation (8), enabling the computation of the temperature at the next time level using current values at the surrounding nodes.

4. Stability Analysis

In the numerical solution of partial differential equations, such as the two-dimensional heat equation, ensuring the stability of the numerical scheme is essential for obtaining reliable and accurate results. Stability determines whether the numerical errors introduced at each time step will remain bounded or amplify uncontrollably as the simulation progresses. Particularly for explicit schemes like the Forward-Time Central-Space (FTCS) method, stability imposes strict conditions on the relationship between time and spatial discretization parameters. Without satisfying these conditions, the computed solution may diverge from the true behavior of the underlying physical system, rendering the results unreliable. In what follows, we conduct a detailed Von Neumann stability analysis of the FTCS scheme applied to the 2D heat equation to derive the necessary condition for numerical stability [14] [15] [19].

Assume the error (or solution) at grid point $\left(i,j\right)$ and time step $n$ is a Fourier mode:

$u_{i,j}^n=G^n{\text{e}}^{\text{i}\left(k_x{x}_i+k_y{y}_j\right)}=G^n{\text{e}}^{\text{i}\left(k_{x}i\text{Δ}x+k_{y}j\text{Δ}y\right)},$

where,

$G$ is the amplification factor;

$k_x,k_y$ are spatial frequencies (wave numbers);

$\text{i}=\sqrt{-1}$ is the imaginary unit.

So each term in FTCS updated formula (8) becomes:

$u_{i,j}^n=G^n{\text{e}}^{\text{i}\left(k_{x}i\text{Δ}x+k_{y}j\text{Δ}y\right)}$

$u_{i,j}^{n+1}=G^{n+1}{\text{e}}^{\text{i}\left(k_{x}i\text{Δ}x+k_{y}j\text{Δ}y\right)}$

$u_{i+1,j}^n=G^n{\text{e}}^{\text{i}\left(k_x\left(i+1\right)\text{Δ}x+k_{y}j\text{Δ}y\right)}=G^n{\text{e}}^{\text{i}{k}_x\text{Δ}x}{\text{e}}^{\text{i}\left(k_{x}i\text{Δ}x+k_{y}j\text{Δ}y\right)}$

$u_{i-1,j}^n=G^n{\text{e}}^{-\text{i}{k}_x\text{Δ}x}{\text{e}}^{\text{i}\left(k_{x}i\text{Δ}x+k_{y}j\text{Δ}y\right)}$

$u_{i,j+1}^n=G^n{\text{e}}^{\text{i}{k}_y\text{Δ}y}{\text{e}}^{\text{i}\left(k_{x}i\text{Δ}x+k_{y}j\text{Δ}y\right)}$

$u_{i,j-1}^n=G^n{\text{e}}^{-\text{i}{k}_y\text{Δ}y}{\text{e}}^{\text{i}\left(k_{x}i\text{Δ}x+k_{y}j\text{Δ}y\right)}.$

Now substitute these values in (8), and factor out the common term $G^n{\text{e}}^{\text{i}\left(k_{x}i\text{Δ}x+k_{y}j\text{Δ}y\right)}$ from the entire right-hand side, we obtain

$\begin{array}{c}{G}^{n+1}{\text{e}}^{\text{i}\left(k_{x}i\text{Δ}x+k_{y}j\text{Δ}y\right)}=G^n{\text{e}}^{\text{i}\left(k_{x}i\text{Δ}x+k_{y}j\text{Δ}y\right)}[1+S_x\left({\text{e}}^{\text{i}{k}_x\text{Δ}x}+{\text{e}}^{-\text{i}{k}_x\text{Δ}x}-2\right)\\ +S_y\left({\text{e}}^{\text{i}{k}_y\text{Δ}y}+{\text{e}}^{-\text{i}{k}_y\text{Δ}y}-2\right)]\end{array}$

Using the identity ${\text{e}}^{\text{i}\theta }+{\text{e}}^{-\text{i}\theta }=2\text{cos}\left(\theta \right)$ , we simplify:

$\begin{array}{l}{G}^{n+1}{\text{e}}^{\text{i}\left(k_{x}i\text{Δ}x+k_{y}j\text{Δ}y\right)}\\ =G^n{\text{e}}^{\text{i}\left(k_{x}i\text{Δ}x+k_{y}j\text{Δ}y\right)}\left[1+2S_x\left(\text{cos}\left(k_x\text{Δ}x\right)-1\right)+2S_y\left(\text{cos}\left(k_y\text{Δ}y\right)-1\right)\right]\end{array}$

$\Rightarrow G=1+2S_x\left(\cos \left(k_x\Delta x\right)-1\right)+2S_y\left(\cos \left(k_y\Delta y\right)-1\right).$

Or, more compactly using the identity $\cos \left(\theta \right)-1=-2{\sin }^2\left(\frac{\theta }{2}\right)$ , we get

$G=1-4S_x{\sin }^2\left(\frac{k_x\Delta x}{2}\right)-4S_y{\sin }^2\left(\frac{k_y\Delta y}{2}\right).$

The scheme is stable if $\left|G\right|\le 1$ for all wave numbers $k_x,k_y$ . Since ${\text{sin}}^2(\cdot )\in \left[0,1\right]$ , the maximum dissipation occurs when

${\sin }^2\left(\frac{k_x\Delta x}{2}\right)=1,{\sin }^2\left(\frac{k_y\Delta y}{2}\right)=1.$

Then the worst-case amplification factor is,

$G_{\text{min}}=1-4S_x-4S_y.$

To ensure $\left|G\right|\le 1$ , the condition is

$-1\le 1-4S_x-4S_y\le 1.$

The most restrictive bound is given by

$\begin{array}{l}1-4S_x-4S_y\ge -1\\ \Rightarrow 4S_x+4S_y\le 2\\ \Rightarrow S_x+S_y\le \frac{1}{2},\end{array}$ (9)

The equality case $S_x+S_y=\frac{1}{2}$ results in marginal stability, where the amplification factor $G=-1$ in the worst-case scenario. This may result in non-growing but oscillatory behavior, which is why it is treated cautiously in practice. Equation (9) is known as the CFL (Courant-Friedrichs-Lewy) condition for stability of the 2D FTCS scheme [20]. This condition ensures the numerical scheme does not produce oscillations or divergence.

To illustrate the behavior of the FTCS scheme under various stability conditions, we present six numerical solution plots for the two-dimensional heat equation using different values of the stability parameter $S_x+S_y$ . Figure 2 and Figure 3 demonstrate solutions within or at the edge of the theoretical stability limit, while Figure 4 shows the unstable behavior when the CFL condition is violated.

Figure 2(a) and Figure 2(b) display stable solutions for $S_x+S_y=0.4$ , simulated over the spatial domains $x\in \left[0,\pi \right]$ and $x\in \left[0,2\pi \right]$ , respectively. Similarly, Figure 3(a) and Figure 3(b) show solutions for $S_x+S_y=0.5$ , which lies exactly at the stability threshold. In all four cases, the numerical solutions remain smooth and bounded over time, indicating stable evolution of heat diffusion.

In contrast, Figure 4(a) and Figure 4(b) illustrate the unstable nature of the FTCS scheme for $S_x+S_y=0.6$ , which exceeds the stability bound. The numerical solutions exhibit unbounded growth and spurious oscillations, especially near the center of the domain, confirming the breakdown of the method under instability.

Note that, the stability parameter $S_x+S_y=0.4$ is computed using $\text{Δ}x=\text{Δ}y=0.05$ , $\text{Δ}t=0.0005$ , and $\alpha =1$ , which yields $S_x=S_y=0.2$ . Similarly, $S_x+S_y=0.5$ results from $\text{Δ}t=0.000625$ , giving $S_x=S_y=0.25$ . On the other hand, $S_x+S_y=0.6$ corresponds to $\text{Δ}t=0.00075$ , yielding $S_x=S_y=0.3$

Figure 2. Stable FTCS solutions of the 2D heat equation for $S_x+S_y=0.4$ .

Figure 3. FTCS solutions at the critical stability threshold $S_x+S_y=0.5$ .

Figure 4. Unstable FTCS solutions of the 2D heat equation for $S_x+S_y=0.6$ .

5. Error Analysis and Results Discussion

To evaluate the accuracy and reliability of the Forward Time Central Space (FTCS) method for solving the two-dimensional heat equation, we perform a quantitative and qualitative error analysis by comparing the numerical solution with the exact analytical solution at a fixed final time $T$ . We chose the final time $T=1$ to ensure the solution has evolved meaningfully while maintaining numerical stability. This duration is long enough to observe the diffusion process without violating the CFL condition [13] [21]. The governing equation is discretized using a uniform grid of size 21 × 21, with spatial step sizes $\text{Δ}x=\text{Δ}y=0.05$ , and time step $\text{Δ}t=0.0005$ . The diffusion coefficient $\alpha$ is set to unity.

The key stability parameter for the FTCS method in two dimensions is the sum $S_x+S_y$ , where

$S_x=\frac{\alpha \text{Δ}t}{\text{Δ}{x}^2},S_y=\frac{\alpha \text{Δ}t}{\text{Δ}{y}^2}.$

In our setup, $\text{Δ}x=\text{Δ}y=0.05$ , $\text{Δ}t=0.0005$ , and $\alpha =1$ , giving $S_x=S_y=0.2$ , yielding $S_x+S_y=0.4$ , which satisfies the theoretical stability criterion (9). This ensures that the scheme remains numerically stable throughout the simulation.

We use the exact solution,

$u_{\text{exact}}\left(x,y,T\right)=\text{sin}\left(\pi x\right)\text{sin}\left(\pi y\right)\text{exp}\left(-2{\pi }^2\alpha T\right)$ (10)

and the sinusoidal initial condition,

$u\left(x,y,0\right)=\text{sin}\left(\pi x\right)\text{sin}\left(\pi y\right)$ (11)

for comparison. The condition (11) is chosen because it satisfies the homogeneous Dirichlet boundary conditions (3 - 4) and corresponds to the exact analytical solution (10) of the 2D heat Equation (1). This facilitates validation of the numerical method by allowing precise error comparison. The error is assessed both visually and through quantitative norms. The maximum norm error is given by

${\Vert e\Vert }_{\infty }=\underset{i,j}{\text{max}}\left|u_{i,j}^{\text{FTCS}}-u_{i,j}^{\text{exact}}\right|\approx 1.48472\times {10}^{-10},$

indicating the largest pointwise deviation across the computational domain.

The discrete $L^2$ -norm error is computed as

${\Vert e\Vert }_2=\sqrt{{\displaystyle \sum }_{i,j}{\left(u_{i,j}^{\text{FTCS}}-u_{i,j}^{\text{exact}}\right)}^2\text{Δ}x\text{Δ}y}\approx 7.42362\times {10}^{-11},$

providing a measure of the global error over the entire domain.

The numerical solution obtained using the FTCS scheme (solid blue line) is plotted against the exact analytical solution (dashed red line) in Figure 5(a). Both solutions are evaluated along the horizontal cross-section $y=0.5$ . The two curves overlap almost perfectly, demonstrating excellent agreement between the FTCS approximation and the true solution. This alignment confirms that the spatial and temporal discretizations are sufficiently fine to resolve the solution accurately.

The second figure of Figure 5 (i.e., Figure 5(b)) plots the absolute pointwise error $\left|u_{\text{FTCS}}-u_{\text{exact}}\right|$ along the same line $y=0.5$ . The error curve exhibits a smooth and symmetric profile, peaking near the center of the domain but remaining under 1.5 × 10⁻¹⁰, which is consistent with round-off level errors given the double-precision arithmetic typically used in MATLAB.

Overall, the error analysis confirms that the FTCS method is both stable and highly accurate under the given discretization and time step. The small magnitude of the maximum and L² errors suggest that the numerical solution is an excellent approximation to the analytical solution, validating the correctness and precision of the implementation.

Figure 5. Comparison of FTCS and Exact Solutions (left) with Corresponding Error (right) at $y=0.5$ , $T=1$ .

6. Conclusions

In this study, we applied the Forward Time Centered Space (FTCS) explicit finite difference scheme to numerically solve the two-dimensional heat conduction equation. Starting from the model formulation, we derived the FTCS update formula through systematic spatial and temporal discretization. A detailed von Neumann stability analysis was performed to obtain the CFL condition, which ensures the boundedness of the numerical solution over time.

The numerical experiments demonstrated that the FTCS method, when applied under stable conditions, yields excellent agreement with the analytical solution. The computed maximum norm and L²-norm errors were both of the order of 10⁻¹⁰, indicating the accuracy and reliability of the method. Graphical comparisons further confirmed the validity of the implementation and the method’s capability to capture the correct behavior of the heat distribution.

A key limitation of FTCS is its conditional stability, requiring very small time steps for fine spatial grids (due to $\text{Δ}t\propto {\left(\text{Δ}x\right)}^2$ ). This leads to high computational cost and makes the scheme unsuitable for long-time simulations or large-scale domains unless adaptive or implicit methods are used. Nevertheless, for problems where the CFL condition can be satisfied, the FTCS method remains a powerful and educationally valuable tool for approximating solutions to the heat equation.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Mojumder, M.S.H., Haque, M.N. and Alam, M.J. (2023) Efficient Finite Difference Methods for the Numerical Analysis of One-Dimensional Heat Equation. Journal of Applied Mathematics and Physics, 11, 3099-3123. https://doi.org/10.4236/jamp.2023.1110204
[2]	Recktenwald, G.W. (2004) Finite Difference Approximations to the Heat Equation. Portland State University.
[3]	Kinkino Meyu, G. and Aliyi Koriche, K. (2021) Analytical Solution vs. Numerical Solution of Heat Equation Flow through Rod of Length 8 Units in One Dimension. International Journal of Applied Mathematics and Theoretical Physics, 7, 53-61. https://doi.org/10.11648/j.ijamtp.20210702.12
[4]	Loskor, W.Z. and Sarkar, R. (2022) A Numerical Solution of Heat Equation for Several Thermal Diffusivity Using Finite Difference Scheme with Stability Conditions. Journal of Applied Mathematics and Physics, 10, 449-465. https://doi.org/10.4236/jamp.2022.102034
[5]	Kabir, M.H., Afroz, A., Hossain, M.M. and Gani, M.O. (2010) A Finite Difference Method for One Dimensional Heat Equation. Journal of Scientific and Engineering Research, 1, 54-59.
[6]	Tawade, S.V. (2022) Comparative Analysis of Two-Dimensional Steady State Heat Flow Problem of a Rectangular Plate by Analytical and Numerical Approach. International Journal of Advanced Research in Science, Communication and Technology, 2, 405-412. https://doi.org/10.48175/ijarsct-2564
[7]	Bernatz, R.A. (2010) Fourier Series and Numerical Methods for Partial Differential Equations. Wiley. https://doi.org/10.1002/9780470651384
[8]	Makhtoumi, M. (2010) Numerical of Solutions Heat Diffusion Equation over One Dimensional Rod Region. International Journal of Advanced Science and Technology, 7, 10-13.
[9]	Gorguis, A. and Benny Chan, W.K. (2008) Heat Equation and Its Comparative Solutions. Computers & Mathematics with Applications, 55, 2973-2980. https://doi.org/10.1016/j.camwa.2007.11.028
[10]	Adomian, G. (1988) A Review of the Decomposition Method in Applied Mathematics. Journal of Mathematical Analysis and Applications, 135, 501-544. https://doi.org/10.1016/0022-247x(88)90170-9
[11]	Adomian, G., Rach, R. and Elrod, M. (1989) On the Solution of Partial Differential Equations with Specified Boundary Conditions. Journal of Mathematical Analysis and Applications, 140, 569-581. https://doi.org/10.1016/0022-247x(89)90084-x
[12]	Mastoi, S., Ganie, A.H., Saeed, A.M., Ali, U., Rajput, U.A. and Mior Othman, W.A. (2022) Numerical Solution for Two-Dimensional Partial Differential Equations Using SM’s Method. Open Physics, 20, 142-154. https://doi.org/10.1515/phys-2022-0015
[13]	Gonçalves de Brito dos Santos, V. and Gomes dos Anjos, P.T. (2022) Finite Difference Method Applied in Two-Dimensional Heat Conduction Problem in the Permanent Regime in Rectangular Coordinates. Advances in Pure Mathematics, 12, 505-518. https://doi.org/10.4236/apm.2022.129038
[14]	Li, N., Steiner, J. and Tang, S. (1994) Convergence and Stability Analysis of an Explicit Finite Difference Method for 2-Dimensional Reaction-Diffusion Equations. The Journal of the Australian Mathematical Society. Series B. Applied Mathematics, 36, 234-241. https://doi.org/10.1017/s0334270000010377
[15]	Rajput, M.N., Shaikh, A.A. and Kamboh, S.A. (2020) Computational Analysis of the Stability of 2D Heat Equation on Elliptical Domain Using Finite Difference Method. Asian Research Journal of Mathematics, 16, 8-19. https://doi.org/10.9734/arjom/2020/v16i330177
[16]	Haque, M.N., Nasu, N.J., Al Mahbub, M.A. and Mohebujjaman, M. (2024) Two-Grid Stabilized Lowest Equal-Order Finite Element Method for the Dual-Permeability-Stokes Fluid Flow Model. Journal of Scientific Computing, 102, Article No. 1. https://doi.org/10.1007/s10915-024-02723-x
[17]	Mahardika, D.P. and Haryani, F.F. (2020) Numerical Analysis of One Dimensional Heat Transfer on Varying Metal. Journal of Physics: Conference Series, 1511, Article 012049. https://doi.org/10.1088/1742-6596/1511/1/012049
[18]	Firmansah, E. and Rosyid, M.F. (2019) Study of Heat Equations with Boundary Differential Equations. Journal of Physics: Conference Series, 1170, Article 012015. https://doi.org/10.1088/1742-6596/1170/1/012015
[19]	Ali, U., Abdullah, F.A. and Ismail, A.I. (2009) Crank-Nicolson Finite Difference Method for Two-Dimensional Fractional Sub-Diffusion Equation. Journal of Interpolation and Approximation in Scientific Computing, 1, 71-77.
[20]	LeVeque, R.J. (2007) Finite Difference Methods for Ordinary and Partial Differential Equations. Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898717839
[21]	Stynes, M., O’Riordan, E. and Gracia, J.L. (2017) Error Analysis of a Finite Difference Method on Graded Meshes for a Time-Fractional Diffusion Equation. SIAM Journal on Numerical Analysis, 55, 1057-1079. https://doi.org/10.1137/16m1082329