Numerical solution of diffusion equation using a method of lines and generalized finite differences

Latest revision as of 08:34, 23 June 2022

Abstract

One of the greatest challenges in the area of applied mathematics continues to be the design of numerical methods capable of approximating the solution of partial differential equations quickly and accurately. One of the most important equations, due to the hydraulic and transport applications it has, and the large number of difficulties that it usually presents when solving it numerically is the Diffusion Equation.

In the present work, a Method of Lines applied to the numerical solution of the said equation in irregular regions is presented using a scheme of Generalized Finite Differences. The second-order finite difference method uses a central node and 8 neighbor points in order to address the spatial approximation. A series of tests and numerical results are presented, which show the accuracy of the proposed method.

Keywords: Generalized finite differences, method of lines, diffusion equation, irregular regions

1. Introduction

When a substance is been transported with a speed equal to zero, i.e. no flow moves it, then, according to the transport equation [1], the initial profile of the substance should remain the same as time goes by, maintaining the substance concentration at each point in the fluid.

However, in practice, this does not happen due to molecular diffusion; this is due to the molecules being in constant motion, causing collisions and rebounds in different directions. The molecules will tend to separate, or diffuse, within the flow. In general, this movement moves from areas with a higher density of molecules to areas with a lower density; this is called diffusion and can be described by a well-known Partial Differential Equation:

(1)

where ${\textstyle \nu }$ represents the diffusion coefficient, which tells how easy a substance diffuses in a medium.

In the past years, many people have worked in different methods to achieve good numerical solutions to this equation, involving a great variety of modern and classical techniques. Nevertheless, even though a large number of numerical methods have been proposed to solve it, a large number of these have a rather limited application to real-world scenarios since they are designed for regular regions.

This is due to the fact that the spatial discretization of the diffusion equation imposes bounds on the time step size in order to achieve numerical stability. However, the use of semidiscretization schemes allows for overcoming many stability issues by using well-known and widely used one-step methods for ordinary initial value problems. For example, in Manshoor et al. [2], a Method of Lines, involving solutions with ${\textstyle a}$ Runge-Kutta method, is presented along with its stability analysis; the results presented show that it is possible to use this kind of method to compute numerical solutions of the equation, yet, the regions where the method is tested are unidimensional regions. On the side of the Generalized Finite Difference methods, Ureña et al. [3] present a scheme to achieve numerical solutions using this method; the results presented in this work show that it is possible to solve different Partial Differential Equations using finite differences over non-regular clouds of points; this work presents that the solutions obtained with this method satisfactorily match with the exact solution of the proposed tests. Some other authors, like Li et al. [4,5,6,7,8,9], and Wang et al. [10] have addressed a large number of applications. Nevertheless, the computational cost of the method as proposed in the aforementioned papers is high, so for each node of the cloud it is required to take up to twenty-six support nodes, even in 2D problems, which increases the computational cost of the implementation.

On the other hand, some variations of the presented generalized finite difference method for several transport equations, which produce satisfactory numerical solutions using a low-cost implementation for the spatial discretization were presented in [11,12,13,14,15]. Even though, the use of several straightforward time integration schemes remained an important issue to take into account.

For the case of interest of this work, it is important to obtain an approximation in generalized finite differences, for the spatial part, to the solution of the problem

where ${\textstyle \Omega }$ is a simply connected planar domain, and its boundary, and ${\textstyle \partial \Omega }$ is a positively oriented Jordan polygon, as the domain shown in Figure 1.

Figure 1. Example of an

domain

On the other hand, for the temporal discretization, a Method of Lines (MOL) is proposed [16]. The basic idea of the MOL is to solve a time-dependent Partial Differential Equation (PDE) by discretizing the spatial derivatives and then, integrating the semi-discretized problem as an Ordinary Differential Equations (ODE) system.

2. Proposed scheme

In order to apply a MOL for the case of the diffusion equation

it is possible to discretize the spatial derivatives applying a generalized finite differences method, for that it is convenient to considerate the approximation to the second or-er linear operator

(2)

where ${\textstyle A}$ , ${\textstyle B}$ , ${\textstyle C}$ , ${\textstyle D}$ , ${\textstyle E}$ , and ${\textstyle F}$ are given functions. Within an arbitrary distribution of nodes, like the one presented in Figure 2, it is possible to approximate its value at a node ${\textstyle p_{0}=(x_{0},y_{0})}$ using values of ${\textstyle u}$ at some neighbor nodes ${\textstyle p_{i}=(x_{i},y_{i})}$ , ${\textstyle i~=~1,2,\dots ,q}$ [12]. For this work, a finite difference scheme is applied at the node ${\textstyle p_{0}}$ , which can be written as a linear combination as

(3)

where ${\textstyle \Gamma _{0},\Gamma _{1},\cdots ,\Gamma _{q}}$ are adequate weighs.

Figure 2. Arbitrary distribution of

and its neighbors

According to Strikwerda [17] and Thomas [18], a finite difference scheme ${\textstyle L_{0}}$ is consistent with the linear operator ${\textstyle L}$ if the local truncation error ${\textstyle \tau }$ satisfies that

(4)

as ${\textstyle p_{1}}$ , ${\textstyle p_{2}}$ , ${\textstyle \dots }$ , ${\textstyle p_{q}\to p_{0}}$ .

Using the six first terms of Taylor's series expansion, up to second order, of the consistency condition (4), it is possible to obtain the system

\left({\begin{array}{cccc}1&1&...&1\\0&\Delta x_{1}&...&\Delta x_{q}\\0&\Delta y_{1}&...&\Delta y_{q}\\0&(\Delta x_{1})^{2}&...&(\Delta x_{q})^{2}\\0&\Delta x_{1}\Delta y_{1}&...&\Delta x_{q}\Delta y_{q}\\0&(\Delta y_{1})^{2}&...&(\Delta y_{q})^{2}\\\end{array}}\right)\left({\begin{array}{c}\Gamma _{0}\\\Gamma _{1}\\\Gamma _{2}\\\cdot \\\cdot \\\cdot \\\Gamma _{q}\\\end{array}}\right)=\left({\begin{array}{c}F(p_{0})\\D(p_{0})\\E(p_{0})\\2A(p_{0})\\B(p_{0})\\2C(p_{0})\\\end{array}}\right),

(5)

where ${\textstyle \Delta x_{i}=x_{i}-x_{0}}$ and ${\textstyle \Delta y_{i}=y_{i}-y_{0}}$ . In order o solve this linear system, it is possible to separate the first equation of the system (5)

(6)

and then, the problem defined by

\left({\begin{array}{ccc}\Delta x_{1}&...&\Delta x_{q}\\\Delta y_{1}&...&\Delta y_{q}\\(\Delta x_{1})^{2}&...&(\Delta x_{q})^{2}\\\Delta x_{1}\Delta y_{1}&...&\Delta x_{q}\Delta y_{q}\\(\Delta y_{1})^{2}&...&(\Delta y_{q})^{2}\\\end{array}}\right)\left({\begin{array}{c}\Gamma _{1}\\\Gamma _{2}\\\cdot \\\cdot \\\cdot \\\Gamma _{q}\\\end{array}}\right)=\left({\begin{array}{c}D(p_{0})\\E(p_{0})\\2A(p_{0})\\B(p_{0})\\2C(p_{0})\\\end{array}}\right),

(7)

can be solved using the reduced Cholesky factorization of its normal equations, as in Tinoco-Guerrero et al. [19], namely

where

M=\left({\begin{array}{ccc}\Delta x_{1}&...&\Delta x_{q}\\\Delta y_{1}&...&\Delta y_{q}\\(\Delta x_{1})^{2}&...&(\Delta x_{q})^{2}\\\Delta x_{1}\Delta y_{1}&...&\Delta x_{q}\Delta y_{q}\\(\Delta y_{1})^{2}&...&(\Delta y_{q})^{2}\\\end{array}}\right),\quad \Gamma =\left({\begin{array}{c}\Gamma _{1}\\\Gamma _{2}\\\cdot \\\cdot \\\cdot \\\Gamma _{q}\\\end{array}}\right),\quad \beta =\left({\begin{array}{c}D(p_{0})\\E(p_{0})\\2A(p_{0})\\B(p_{0})\\2C(p_{0})\\\end{array}}\right).

The value of ${\textstyle \Gamma _{0}}$ is then obtained from Eq.(6) assuming, for the case of the diffusion equation, ${\textstyle F(p_{0})=0}$ .

Now, the scheme defined by Eq.(7) can be used to approximate the linear operator

taking ${\textstyle A=C=\nu }$ , and ${\textstyle B=D=E=0}$ . The resulting ${\textstyle \Gamma _{i}}$ coefficients, define the Finite-Difference Scheme

(8)

for diffusion equation, where ${\textstyle u_{i}(t)}$ is the approximation to the solution in the point ${\textstyle p_{i}=(x_{i},y_{i})}$ , and ${\textstyle p_{j}}$ are the corresponding neighbor nodes of ${\textstyle p_{i}}$ .

An important issue to be taken into account is the number of neighbors, ${\textstyle q}$ , to use in the scheme. In this paper, 8 neighbor points were taken into account following the stencil shown in Figure 3.

Figure 3. Centered stencil used in the scheme

Once a discretization of the spatial operator is obtained, following the idea presented in Tinoco-Guerrero et al. [1], the semi-discretized PDE can be rewritten as a linear ODE system in time,

where ${\textstyle m}$ is the total amount of grid points.

There exist several ways to approximate the solution of this ODE system, for example, it is possible to use a forward Euler method as

nevertheless, this kind of implementation has proven to be conditionally stable and, in some cases, the stability conditions are difficult to accomplish.

Another way to solve this system is the Method of Lines, where each ${\textstyle u_{i}(t)}$ is solved for a fixed grid node ${\textstyle p_{i}=(x_{i},y_{i})}$ , i.e. it is solved by “lines”. As we are now dealing with a set of ODEs, it is possible to use a Runge-Kutta method, to solve the ODE on each line in order to obtain stable and accurate results.

In the present work it is proposed to use Runge-Kutta [20] of 2nd, 3rd, and 4th order, as follows.

Second-order Runge-Kutta method

where

Third-order Runge-Kutta method

where

Fourth-order Runge-Kutta method

where

in all the cases, ${\textstyle k}$ represents the time level.

3. Numerical tests

To show the performance of the proposed scheme, the problem of obtaining a numerical solution to the diffusion equation in four different regions was proposed. The first region, denoted as A, corresponds to the unitary square, for comparison reasons, and the other three regions, denoted as B, C, and D, are non-rectangular planar domains. All the regions were scaled to fit on ${\textstyle [0,1]\times [0,1]}$ , and meshed with ${\textstyle 21\times {21}}$ nodes, following a variational procedure implemented in UNAMalla [21], then they were subdivided to obtain meshes with ${\textstyle 41\times {41}}$ nodes. The meshes for with ${\textstyle 21\times {21}}$ nodes for each region can be seen in Figure 4.



Figure 4. Meshes of the test regions with $21\times {21}$ nodes

For all the regions, two different tests were performed:

Test 1. Following the same idea as in Tinoco-Ruiz [22], where a well-known diffusion problem is presented, the initial and boundary conditions were taken from the closed-form solution:

Test 2. For comparison purposes, one of the tests presented in Sánchez et al. [23], was also taken into account; in this case, the initial and boundary conditions were taken from the closed-form solution:

The time interval ${\textstyle [0,1]}$ was subdivided with different discretizations, the number of time steps was chosen to satisfy the classical Courant condition [24]:

where ${\textstyle \Delta x}$ and ${\textstyle \Delta y}$ are the spatial steps on each grid, taken from the meshes with ${\textstyle 41\times {41}}$ nodes, and ${\textstyle \nu =0.2}$ .

The norm of the quadratic error, at a ${\textstyle k-}$ the time level can be computed as

where ${\textstyle U_{i,j}^{k}}$ and ${\textstyle u_{i,j}^{k}}$ are the approximated and theoretical solutions, respectively, at the ${\textstyle i,j}$ -th grid node, and ${\textstyle A_{i,j}}$ is the area of the polygon defined by the points ${\textstyle P_{i+1,j}}$ , ${\textstyle P_{i,j+1}}$ , ${\textstyle P_{i-1,j}}$ and ${\textstyle P_{i,j-1}}$ .

Figure 5 presents a comparison of the numerical results obtained with the proposed scheme using the second-order Runge-Kutta approximation and the exact solutions for the test region A at time ${\textstyle t=0.66}$ s. The approximation is presented on the left, and the exact solution on the right.

Figure 5. Comparison between numerical results (left) and exact solution (right) for test region A. Test 1 (Top) and Test 2 (Bottom)

Similarly, Figures 6, 7, and 8 present the respective comparisons for the regions B, C, and D.

Figure 6. Comparison between numerical results (left) and exact solution (right) for test region B. Test 1 (Top) and Test 2 (Bottom)

Figure 7. Comparison between numerical results (left) and exact solution (right) for test region C. Test 1 (Top) and Test 2 (Bottom)

Figure 8. Comparison between numerical results (left) and exact solution (right) for test region D. Test 1 (Top) and Test 2 (Bottom)

Figures 9 to 12 present the maximum value of the error ${\textstyle \|e^{k}\|_{2}}$ computed for all the regions. For each figure, both tests were performed in the regions meshed with ${\textstyle 21\times {21}}$ nodes and ${\textstyle 41\times {41}}$ nodes.


Figure 9. Maximum $\\|e^{k}\\|_{2}$ computed for test region A with $21\times {21}$ nodes (Left) and $41\times {41}$ nodes (Right)


Figure 10. Maximum $\\|e^{k}\\|_{2}$ computed for test region B with $21\times {21}$ nodes (Left) and $41\times {41}$ nodes (Right)


Figure 11. Maximum $\\|e^{k}\\|_{2}$ computed for test region C with $21\times {21}$ nodes (Left) and $41\times {41}$ nodes (Right)


Figure 12. Maximum $\\|e^{k}\\|_{2}$ computed for test region D with $21\times {21}$ nodes (Left) and $41\times {41}$ nodes (Right)

4. Conclusions

The numerical results show that the proposed method of Lines applied to the diffusion equation produces satisfactory numerical solutions. In the tests carried out, no spurious oscillations or instabilities were perceived. In addition, the results show that it is not necessary to “work more” by using higher Runge-Kutta methods; this scheme accomplishes better results when second-order Runge-Kutta is used to solve the Ordinary Differential Equations system.

It is possible to appreciate that with the presented discretizations, even with the fourth-order Runge-Kutta method, acceptable numerical results can be obtained. Additionally, the numerical results shown in Figures 9 to 12 show that a refinement of the spatial mesh can improve the approximations carried out with the method; as expected in these kinds of methods. It is worth mentioning that, for all the tests, a minimum number of time steps was used, taking into account the classical theory for the diffusion equation, which makes the computational cost of this method low, since it does not require making very small temporary discretizations, as in other cases.

Furthermore, the proposed scheme can produce stable results for non-standard initial conditions in highly irregular regions. For example, in the following videos

solutions of the diffusion equation are presented on domains that are geometrical approximations of real geographical locations, where the boundary conditions are fixed as ${\textstyle 0}$ and the initial condition is stated as:

It is possible to see that even for these conditions in these regions, the scheme produces stable results that show the expected behavior.

An important remark is that, even when in this work structured convex grids were used, the development of the method doesn't take into account a particular data structure, i.e. this method can be used not only on structured meshes but also as a meshless method, that would be developed as future work.

Acknowledgments

We want to thank AULA CIMNE-Morelia and CIC-UMSNH for the financial support for this work. Thanks to Universidad Vasco de Quiroga for lending the facilities to perform the final part of the work done for the paper.

Special thanks to the ``Laboratorio de Matemáticas Aplicadas of the ``Facultad de Ciencias Físico-Matemáticas at ``Universidad Michoacana de San Nicolás de Hidalgo for lending their computational infrastructure.

The authors would also like to thank the reviewers for their valuable feedback and suggestions to improve the present work.

References

[1] Tinoco-Guerrero G., Domínguez-Mota F.J., Tinoco-Ruiz J.G., Lucas-Martínez J.S. Aproximación de la ecuación de advección en regiones irregulares utilizando un método de líneas y diferencias finitas generalizadas. Revista Mexicana de Métodos Numéricos, Volume 2, 3, 2028.

[2] Manshoor B., Salleh H., Khalid A., Sayed Abdelaal M.A. Method of lines and Runge-Kutta method in solving partial differential equation for eeat equation. Journal of Complex Flow, 3(1):21-25, 2021.

[3] Ureña Prieto F., Benito Muñoz J.J., Gavete Corvinos L. Application of the generalized finite difference method to solve the advection-diffusion equation. Journal of Computational and Applied Mathematics, 25:1849-1855, 2011.

[4] Li P.-W., Fan C.-M., Grabski J.K. A meshless generalized finite difference method for solving shallow water equations with the flux limiter technique. Engineering Analysis with Boundary Elements, 131:159–173, 2021.

[5] Zhang T., Lin Z.-H., Huang G.-Y., Fan C.-M., Li P.-W. Solving Boussinesq equations with a meshless finite difference method. Ocean Engineering, 198:106957, 2020.

[6] Fu Z.-J., Tang Z.-C., Zhao H.-T., Li P.-W., Rabczuk T. Numerical solutions of the coupled unsteady nonlinear convection-diffusion equations based on generalized finite difference method. The European Physical Journal Plus, 134(6):1–20, 2019.

[7] Zhang T., Ren Y.-F., Yang Z.-Q., Fan C.-M., Li P.-W. Application of generalized finite difference method to propagation of nonlinear water waves in numerical wave flume. Ocean Engineering, 123:278–290, 2016.

[8] Li P.-W., Fan C.-M. Generalized finite difference method for two-dimensional shallow water equations. Engineering Analysis with Boundary Elements, 80:58–71, 2017.

[9] Fu Z.-J., Zhang J., Li P.-W., Zheng J.-H. A semi-Lagrangian meshless framework for numerical solutions of two-dimensional sloshing phenomenon. Engineering Analysis with Boundary Elements, 112:58–67, 2020.

[10] Wang Q., Kim P., Qu W. A hybrid localized meshless method for the solution of transient groundwater flow in two dimensions. Mathematics, 10(3):1-14, 2022.

[11] Domínguez-Mota F.J., Mendoza-Armenta S., Tinoco-Ruiz J.G. Finite difference schemes satisfying an optimality condition. IMACS Series in Computational and Applied Mathematics: MASCOT, 2011.

[12] Tinoco-Guerrero G., Domínguez-Mota F.J., Tinoco-Ruiz J.G. Numerical solution of advection equations using structured grids on plane irregular regions with a finite differences scheme. Boletín de la Sociedad Mexicana de Computación Científica y Sus Aplicaciones, 1(3):29-32, 2015.

[13] Tinoco-Ruiz J.G., Domínguez-Mota F.J., Tinoco-Guerrero G., Fernández-Valdés P.M., Mendoza-Armenta S. Numerical solution of differential equations in irregular plane regions using quality structured convex grids. International Journal of Modeling, Simulation and Scientific Computing, 4(2):1340005, 2013.

[14] Domínguez-Mota F.J., Lucas-Martínez J.S., Tinoco-Guerrero G. A generalized finite-differences scheme used in modeling of a direct and an inverse problem of advection-diffusion. International Journal of Applied Mathematics, 33(4):599-608, 2020.

[15] Chávez-Negrete C., Domínguez-Mota F.J., Santana-Quinteros D. Numerical solution of Richards' equation of water flow by generalized finite differences. Computers and Geotechnics, 101:168–175, 2018.

[16] Schiesser W. The numerical method of lines. San Diego: Academic Press, 1991.

[17] Strikwerda J.C. Finite difference schemes and partial differential equations. Society for Industrial and Applied Mathematics, pp. 427, 2004.

[18] Thomas J.W. Numerical partial differential equations: finite difference methods. Texts in Applied Mathematics (TAM), Volume 22, Springer, 1995.

[19] Tinoco-Guerrero G., Domínguez-Mota F.J., Tinoco-Ruiz J.G. A study of the stability for a generalized finite-difference scheme applied to the advection-diffusion equation. Mathematics and Computer in Simulations, 176:301-311, 2020.

[20] Butcher J.C. Conference on the numerical solution of differential equations. Lecture Notes in Mathematics (LNM), Springer, 109:133-139, 1969.

[21] UNAMalla. An automatic package for numerical grid generation, 2011.

[22] Tinoco-Ruiz J.G., Domínguez-Mota F.J., Tinoco-Guerrero G., Lucas-Martínez J.S. Study of the stability of a generalized finite difference scheme applied to the diffusion equation in irregular 2D space regions using convex grids. IMACS Series in Computational and Applied Mathematics: MASCOT2018, 2018.

[23] Sánchez L.M., Ureña F., Benito J.J., Gavete L. Resolución de la ecuación de difusión en 2-D y 3-D utilizando diferencias finitas generalizadas. Consistencia y Estabilidad. XXI Congreso de Ecuaciones Diferenciales y Aplicaciones, Ciudad Real, pp. 1-8, 2009.

[24] Abbott M.B., McCowan A.D., Warren I.R. Accuracy of short-wave numerical models. Journal of Hydraulic Engineering, 110(10):1287-1301, 1984.

Latest revision as of 08:34, 23 June 2022

Abstract

1. Introduction

2. Proposed scheme

Second-order Runge-Kutta method

Third-order Runge-Kutta method

Fourth-order Runge-Kutta method

3. Numerical tests

4. Conclusions

Acknowledgments

References

Document information

Document Score

Share this document

Keywords

@@ Line 4: / Line 4: @@
 '''G.&nbsp;Tinoco-Guerrero Corresponding author: gerardo.tinoco@umich.mx<sup>a,b</sup>, F.&nbsp;J.&nbsp;Domínguez-Mota<sup>a</sup>, J.&nbsp;A.&nbsp;Guzmán-Torres<sup>a</sup>, J.&nbsp;G.&nbsp;Tinoco-Ruiz<sup>a</sup>'''
 -->
 ==Abstract==
 One of the greatest challenges in the area of applied mathematics continues to be the design of numerical methods capable of approximating the solution of partial differential equations quickly and accurately. One of the most important equations, due to the hydraulic and transport applications it has, and the large number of difficulties that it usually presents when solving it numerically is the Diffusion Equation.
-In the present work, a Method of Lines applied to the numerical solution of the said equation in irregular regions is presented using a scheme of Generalized Finite Differences. A series of tests and numerical results are presented, which show the accuracy of the proposed method.
+In the present work, a Method of Lines applied to the numerical solution of the said equation in irregular regions is presented using a scheme of Generalized Finite Differences. The second-order finite difference method uses a central node and 8 neighbor points in order to address the spatial approximation. A series of tests and numerical results are presented, which show the accuracy of the proposed method.
-'''keywords''' Generalized Finite Differences, Method of Lines, Diffusion Equation, Irregular Regions
+'''Keywords''': Generalized finite differences, method of lines, diffusion equation, irregular regions
-==1 Introduction==
+==1. Introduction==
 When a substance is been transported with a speed equal to zero, i.e. no flow moves it, then, according to the transport equation <span id='citeF-1'></span>[[#cite-1|[1]]], the initial profile of the substance should remain the same as time goes by, maintaining the substance concentration at each point in the fluid.
-However, in practice, this does not happen due to molecular diffusion; this is due to the molecules are in constant motion, causing collisions and rebounds in different directions. The molecules will tend to separate, or diffuse, within the flow. In general, this movement moves from areas with a higher density of molecules to areas with a lower density.
+However, in practice, this does not happen due to molecular diffusion; this is due to the molecules being in constant motion, causing collisions and rebounds in different directions. The molecules will tend to separate, or diffuse, within the flow. In general, this movement moves from areas with a higher density of molecules to areas with a lower density; this is called ''diffusion'' and can be described by a well-known Partial Differential Equation:
-Taking into account that the diffusive flux of a substance, <math display="inline">f(u)</math> in any part of a system is proportional to the gradient of its density, then it can be written as:
 <span id="eq-1"></span>
@@ Line 27: / Line 24: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>f(u) = -\nu \nabla u(x,t) = -\nu \frac{\partial }{\partial x} u(x,t)       </math>
+| style="text-align: center;" | <math>\frac{\partial u}{\partial t} = \nu \nabla ^2 u       </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (1)
@@ Line 34: / Line 31: @@
 where <math display="inline">\nu </math> represents the diffusion coefficient, which tells how easy a substance diffuses in a medium.
-Since the change in total mass in a given interval <math display="inline">[x1, x2]</math> is due to the flux at the extreme points, then
+In the past years, many people have worked in different methods to achieve good numerical solutions to this equation, involving a great variety of modern and classical techniques. Nevertheless, even though a large number of numerical methods have been proposed to solve it,  a large number of these have a rather limited application to real-world scenarios since they are designed for regular regions.
-{| class="formulaSCP" style="width: 100%; text-align: left;"
+This is due to the fact that the spatial discretization of the diffusion equation imposes bounds on the time step size in order to achieve numerical stability. However, the use of semidiscretization schemes allows for overcoming many stability issues by using well-known and widely used one-step methods for ordinary initial value problems. For example, in  Manshoor et al. <span id='citeF-2'></span>[[#cite-2|[2]]], a Method of Lines, involving solutions with <math display="inline"> a </math> Runge-Kutta method, is presented along with its stability analysis; the results presented show that it is possible to use this kind of method to compute numerical solutions of the equation, yet, the regions where the method is tested are unidimensional regions. On the side of the Generalized Finite Difference methods, Ureña et al. <span id='citeF-3'></span>[[#cite-3|[3]]] present a scheme to achieve numerical solutions using this method; the results presented in this work show that it is possible to solve different Partial Differential Equations using finite differences over non-regular clouds of points; this work presents that the solutions obtained with this method satisfactorily match with the exact solution of the proposed tests. Some other authors, like  Li et al.  <span id='citeF-4'></span><span id='citeF-5'></span><span id='citeF-6'></span><span id='citeF-7'></span><span id='citeF-8'></span><span id='citeF-9'></span>[[#cite-4|[4,5,6,7,8,9]]], and Wang et al. <span id='citeF-10'></span>[[#cite-10|[10]]] have addressed a large number of applications. Nevertheless, the computational cost of the method as proposed in the aforementioned papers is high, so for each node of the cloud it is required to take up to twenty-six support nodes, even in 2D problems, which increases the computational cost of the implementation.
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>    \frac{d}{dt}\int _{x_1}^{x_2} u(x,t) dx = f(u(x_1,t)) - f(u(x_2,t)).   </math>
-|}
-|}
-Assuming enough smoothness for the functions <math display="inline">u</math> and <math display="inline">f</math>, the expression can be written as
+On the other hand, some variations of the presented generalized finite difference method for several transport equations, which produce satisfactory numerical solutions using a low-cost implementation for the spatial discretization were presented in <span id='citeF-11'></span><span id='citeF-12'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-15'></span>[[#cite-11|[11,12,13,14,15]]]. Even though, the use of several straightforward time integration schemes remained an important issue to take into account.
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>    \int _{x_1}^{x_2}\left[\frac{\partial }{\partial t}u(x,t) + \frac{\partial }{\partial x}f(u(x,t))\right]dx = 0   </math>
-|}
-|}
-Since this integral must be equal to zero for any <math display="inline">x_1</math> and <math display="inline">x_2</math> values, then the integrand must be identically zero. This is,
-<span id="eq-2"></span>
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>\frac{\partial }{\partial t}u(x,t) + \frac{\partial }{\partial x}f(u(x,t)) = 0,       </math>
-|}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
-|}
-which is a formula know as a conservation law.
-Replacing the flow equation ([[#eq-1|1]]) in the conservation law ([[#eq-2|2]]) it is possible to obtain,
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>    \frac{\partial }{\partial t} u(x,t) + \frac{\partial }{\partial x}\left(-\nu \frac{\partial }{\partial } u(x,t)\right)= 0   </math>
-|}
-|}
-which, making some pertinent derivations, takes the form
-<span id="eq-3"></span>
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>\frac{\partial }{\partial t}u(x,t) = \nu \frac{\partial ^2}{\partial x^2}u(x,t),       </math>
-|}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
-|}
-which is known as the diffusion equation in a spatial dimension. This equation is also known as the heat equation, since heat diffuses in a similar way.
-This same idea can be extended <span id='citeF-2'></span>[[#cite-2|[2]]], following a similar process to the previous one, to define the diffusion equation in a plane as
-<span id="eq-4"></span>
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>\frac{\partial }{\partial t}u(x,y,t) = \nu \left(\frac{\partial ^2}{\partial x^2}u(x,y,t) + \frac{\partial ^2}{\partial y^2}u(x,y,t)\right).       </math>
-|}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
-|}
 For the case of interest of this work, it is important to obtain an approximation in generalized finite differences, for the spatial part, to the solution of the problem
@@ Line 121: / Line 48: @@
 | style="text-align: center;" | <math>    u(x,y,0) = g(x,y)                          (x, y) \in \Omega ,</math>
 |-
-| style="text-align: center;" | <math>    u(x,y,t)|_{\partial \Omega } = g(x,y,t)                 (x, y) \in \Omega  t\in [0,T],   </math>
+| style="text-align: center;" | <math>    u(x,y,t)\Big|_{\partial \Omega } = g(x,y,t)                 (x, y) \in \Omega  t\in [0,T],   </math>
 |}
 |}
-Where <math display="inline">\Omega </math> is a simply connected planar domain, and its boundary, <math display="inline">\partial \Omega </math>, is a positively oriented Jordan polygon, as the domain shown in Figure [[#img-1|1]].
+where <math display="inline">\Omega </math> is a simply connected planar domain, and its boundary, and <math display="inline">\partial \Omega </math> is a positively oriented Jordan polygon, as the domain shown in  [[#img-1|Figure 1]].
-In this case, for the temporal discretization, a Method of Lines is proposed <span id='citeF-3'></span>[[#cite-3|[3]]].
 <div id='img-1'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+{| style="text-align: center; border: 1px solid; margin: 1em auto; width: auto;max-width: auto;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig01.png|300x300px|Example of an Ω domain.]]
+|[[Image:Tinoco_Guerrero_et_al_2022b_9158_fig01.png|316px|Example of an Ω domain.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 1:''' Example of an <math>\Omega </math> domain.
+| colspan="1"| '''Figure 1'''. Example of an <math>\Omega </math> domain
 |}
-==2 Proposed Scheme==
-For the case of diffusion equation
+On the other hand, for the temporal discretization, a Method of Lines (MOL) is proposed <span id='citeF-16'></span>[[#cite-16|[16]]]. The basic idea of the MOL is to solve a time-dependent Partial Differential Equation (PDE) by discretizing the spatial derivatives and then, integrating the semi-discretized problem as an Ordinary Differential Equations (ODE) system.
+==2. Proposed scheme==
+In order to apply a MOL for the case of the diffusion equation
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 150: / Line 78: @@
 |}
-it is possible to discretize the spatial derivatives applying a generalized finite differences method, for that it is convenient to considerate the approximation to the second order linear operator
+it is possible to discretize the spatial derivatives applying a generalized finite differences method, for that it is convenient to considerate the approximation to the second or-er linear operator
-<span id="eq-5"></span>
+<span id="eq-2"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 160: / Line 88: @@
 | style="text-align: center;" | <math>Lu = Au_{xx} + Bu_{xy} + Cu_{yy} + Du_{x} + Eu_{y} + Fu,       </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
 |}
-where <math display="inline">A</math>, <math display="inline">B</math>, <math display="inline">C</math>, <math display="inline">D</math>, <math display="inline">E</math>, and <math display="inline">F</math> are given functions. Within an arbitrary distribution of nodes, as the one presented on Figure [[#img-2|2]], it is possible to approximate its value at a node <math display="inline">p_0 = (x_0,y_0)</math> using values of <math display="inline">u</math> at some neighbor nodes <math display="inline">p_i = (x_i, y_i)</math>, <math display="inline">i~=~1,2,\dots ,q</math> <span id='citeF-4'></span>[[#cite-4|[4]]]. For this work, a finite difference scheme is applied at the node <math display="inline">p_0</math>, which can be written as a linear combination as
+where <math display="inline">A</math>, <math display="inline">B</math>, <math display="inline">C</math>, <math display="inline">D</math>, <math display="inline">E</math>, and <math display="inline">F</math> are given functions. Within an arbitrary distribution of nodes, like the one presented in  [[#img-2|Figure 2]], it is possible to approximate its value at a node <math display="inline">p_0 = (x_0,y_0)</math> using values of <math display="inline">u</math> at some neighbor nodes <math display="inline">p_i = (x_i, y_i)</math>, <math display="inline">i~=~1,2,\dots ,q</math> <span id='citeF-12'></span>[[#cite-12|[12]]]. For this work, a finite difference scheme is applied at the node <math display="inline">p_0</math>, which can be written as a linear combination as
-<span id="eq-6"></span>
+<span id="eq-3"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 173: / Line 101: @@
 | style="text-align: center;" | <math>L_0 = \Gamma _0 u(p_0) + \Gamma _1 u(p_1) + ... + \Gamma _q u(p_q) = \sum _{i=0}^{q}\Gamma _i u(p_i).   </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
 |}
-where <math display="inline">\Gamma _0, \Gamma _1, \dots , \Gamma _q</math> are adequate weighs.
+where <math display="inline">\Gamma _0, \Gamma _1, \cdots , \Gamma _q</math> are adequate weighs.
 <div id='img-2'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+{| style="text-align: center; border: 1px solid; margin: 1em auto; width: auto;max-width: auto;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig02.png|300x300px|Arbitrary distribution of p₀ and its neghbors.]]
+|[[Image:Tinoco_Guerrero_et_al_2022b_5181_fig02.png|283px|Arbitrary distribution of p₀ and its neighbors.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 2:''' Arbitrary distribution of <math>p_0</math> and its neghbors.
+| colspan="1" | '''Figure 2'''. Arbitrary distribution of <math>p_0</math> and its neighbors
 |}
-According to <span id='citeF-5'></span><span id='citeF-6'></span>[[#cite-5|[5,6]]], a finite difference scheme <math display="inline">L_0</math> is consistent with the linear operator <math display="inline">L</math> if the local truncation error <math display="inline">\tau </math> satisfies that
-<span id="eq-7"></span>
+According to Strikwerda [[#cite-17|[17]]] and Thomas <span id='citeF-18'></span>[[#cite-18|[18]]], a finite difference scheme <math display="inline">L_0</math> is consistent with the linear operator <math display="inline">L</math> if the local truncation error <math display="inline">\tau </math> satisfies that
+<span id="eq-4"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 196: / Line 125: @@
 | style="text-align: center;" | <math>\tau = \left[Lu\right]_{p_{0}} - [L_{0}u]_{p_{0}} \to 0       </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
 |}
 as <math display="inline">p_{1}</math>, <math display="inline">p_{2}</math>, <math display="inline">\dots </math>, <math display="inline">p_{q} \to p_{0}</math>.
-Using the six first terms of Taylor's series expansion, up to second order, of the consistency condition ([[#eq-7|7]]), it is possible to obtain the system
+Using the six first terms of Taylor's series expansion, up to second order, of the consistency condition ([[#eq-4|4]]), it is possible to obtain the system
-<span id="eq-8"></span>
+<span id="eq-5"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 209: / Line 138: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>\left(     \begin{array}{lclc}1 & 1    & ... & 1\\      0 & \Delta x_1   & ... & \Delta x_q\\      0 & \Delta y_1   & ... & \Delta y_q\\      0 & (\Delta x_1)^2  & ... & (\Delta x_q)^2\\      0 & \Delta x_1\Delta y_1 & ... & \Delta x_q \Delta y_q\\      0 & (\Delta y_1)^2  & ... & (\Delta y_q)^2\\     \end{array}    \right)    \left(     \begin{array}{l}\Gamma _0\\      \Gamma _1\\      \Gamma _2\\      .\\      .\\      .\\      \Gamma _q\\     \end{array}    \right)    =     \left(    \begin{array}{l}F(p_0)\\       D(p_0)\\      E(p_0)\\      2A(p_0)\\      B(p_0)\\      2C(p_0)\\     \end{array}    \right),       </math>
+| style="text-align: center;" | <math>\left(     \begin{array}{cccc}1 & 1    & ... & 1\\      0 & \Delta x_1   & ... & \Delta x_q\\      0 & \Delta y_1   & ... & \Delta y_q\\      0 & (\Delta x_1)^2  & ... & (\Delta x_q)^2\\      0 & \Delta x_1\Delta y_1 & ... & \Delta x_q \Delta y_q\\      0 & (\Delta y_1)^2  & ... & (\Delta y_q)^2\\     \end{array}    \right)    \left(     \begin{array}{c}\Gamma _0\\      \Gamma _1\\      \Gamma _2\\      \cdot\\      \cdot\\      \cdot\\      \Gamma _q\\     \end{array}    \right)    =     \left(    \begin{array}{c}F(p_0)\\       D(p_0)\\      E(p_0)\\      2A(p_0)\\      B(p_0)\\      2C(p_0)\\     \end{array}    \right),       </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
 |}
-where <math display="inline">\Delta x_{i} = x_{i} - x_{0}</math> and <math display="inline">\Delta y_{i} = y_{i} - y_{0}</math>. In order to solve this linear system, it is possible to separate the first equation of the system ([[#eq-8|8]])
+where <math display="inline">\Delta x_{i} = x_{i} - x_{0}</math> and <math display="inline">\Delta y_{i} = y_{i} - y_{0}</math>. In order o solve this linear system, it is possible to separate the first equation of the system ([[#eq-5|5]])
-<span id="eq-9"></span>
+<span id="eq-6"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 224: / Line 153: @@
 | style="text-align: center;" | <math>\sum _{i=0}^{q}\Gamma _{i} - F(p_0) = 0       </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
 |}
 and then, the problem defined by
-<span id="eq-10"></span>
+<span id="eq-7"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 235: / Line 164: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>\left(     \begin{array}{lcl}\Delta x_1   & ... & \Delta x_q\\      \Delta y_1   & ... & \Delta y_q\\      (\Delta x_1)^2  & ... & (\Delta x_q)^2\\      \Delta x_1\Delta y_1 & ... & \Delta x_q \Delta y_q\\      (\Delta y_1)^2  & ... & (\Delta y_q)^2\\     \end{array}    \right)    \left(     \begin{array}{l}\Gamma _1\\      \Gamma _2\\      .\\      .\\      .\\      \Gamma _q\\     \end{array}    \right)    =    \left(     \begin{array}{l}D(p_0)\\      E(p_0)\\      2A(p_0)\\      B(p_0)\\      2C(p_0)\\     \end{array}    \right),       </math>
+| style="text-align: center;" | <math>\left(     \begin{array}{ccc}\Delta x_1   & ... & \Delta x_q\\      \Delta y_1   & ... & \Delta y_q\\      (\Delta x_1)^2  & ... & (\Delta x_q)^2\\      \Delta x_1\Delta y_1 & ... & \Delta x_q \Delta y_q\\      (\Delta y_1)^2  & ... & (\Delta y_q)^2\\     \end{array}    \right)    \left(     \begin{array}{c}\Gamma _1\\      \Gamma _2\\      \cdot\\     \cdot\\      \cdot\\      \Gamma _q\\     \end{array}    \right)    =    \left(     \begin{array}{c}D(p_0)\\      E(p_0)\\      2A(p_0)\\      B(p_0)\\      2C(p_0)\\     \end{array}    \right),       </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
 |}
-can be solved using the reduced Cholesky factorization of its normal equations, as in <span id='citeF-7'></span>[[#cite-7|[7]]], namely
+can be solved using the reduced Cholesky factorization of its normal equations, as in Tinoco-Guerrero et al. <span id='citeF-19'></span>[[#cite-19|[19]]], namely
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 258: / Line 187: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>    M = \left(     \begin{array}{clc}      \Delta x_1   & ... & \Delta x_q\\      \Delta y_1   & ... & \Delta y_q\\      (\Delta x_1)^2  & ... & (\Delta x_q)^2\\      \Delta x_1\Delta y_1 & ... & \Delta x_q \Delta y_q\\      (\Delta y_1)^2  & ... & (\Delta y_q)^2\\     \end{array}    \right),    \quad     \Gamma = \left(     \begin{array}{l}      \Gamma _1\\      \Gamma _2\\      .\\      .\\      .\\      \Gamma _q\\     \end{array}    \right),    \quad     \beta = \left(     \begin{array}{l}      D(p_0)\\      E(p_0)\\      2A(p_0)\\      B(p_0)\\      2C(p_0)\\     \end{array}    \right).   </math>
+| style="text-align: center;" | <math>    M = \left(     \begin{array}{ccc}      \Delta x_1   & ... & \Delta x_q\\      \Delta y_1   & ... & \Delta y_q\\      (\Delta x_1)^2  & ... & (\Delta x_q)^2\\      \Delta x_1\Delta y_1 & ... & \Delta x_q \Delta y_q\\      (\Delta y_1)^2  & ... & (\Delta y_q)^2\\     \end{array}    \right),    \quad     \Gamma = \left(     \begin{array}{c}      \Gamma _1\\      \Gamma _2\\      \cdot\\     \cdot\\    \cdot\\      \Gamma _q\\     \end{array}    \right),    \quad     \beta = \left(     \begin{array}{c}      D(p_0)\\      E(p_0)\\      2A(p_0)\\      B(p_0)\\      2C(p_0)\\     \end{array}    \right).   </math>
 |}
 |}
-The value of <math display="inline">\Gamma _0</math> is then obtained from ([[#eq-9|9]]) assuming, for the case of the diffusion equation, <math display="inline">F(p_0) = 0</math>.
+The value of <math display="inline">\Gamma _0</math> is then obtained from Eq.([[#eq-6|6]]) assuming, for the case of the diffusion equation, <math display="inline">F(p_0) = 0</math>.
-Now, the scheme defined by ([[#eq-10|10]]) can be used to approximate the linear operator
+Now, the scheme defined by Eq.([[#eq-7|7]]) can be used to approximate the linear operator
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 277: / Line 206: @@
 taking <math display="inline">A = C = \nu </math>, and <math display="inline">B = D = E = 0</math>. The resulting <math display="inline">\Gamma _{i}</math> coefficients, define the Finite-Difference Scheme
-<span id="eq-11"></span>
+<span id="eq-8"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 283: / Line 212: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>u_{0}^{k+1} = u_{0}^{k} + \sum _{i=0}^{q}\Gamma _{i}u_{i}^{k},       </math>
+| style="text-align: center;" | <math>\frac{\partial }{\partial t}u_{i}(t) = \sum _{j=0}^{q}\Gamma _{j}u_{j}(t),       </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
 |}
-for diffusion equation, where <math display="inline">k</math> represent the time level and <math display="inline">u_{i}^{k}</math> is the approximation to the solution in the point <math display="inline">p_{i} = (x_{i}, y_{i})</math> at time <math display="inline">k\Delta t</math>.
+for diffusion equation, where <math display="inline">u_{i}(t)</math> is the approximation to the solution in the point <math display="inline">p_{i} = (x_{i}, y_{i})</math>, and <math display="inline">p_j</math> are the corresponding neighbor nodes of <math display="inline">p_i</math>.
-An important issue to be taken into account is the number of neighbors, <math display="inline">q</math>, to use in the scheme. In this paper, 8 neighbor points were taken into account following the stencil shown in Figure [[#img-3|3]].
+An important issue to be taken into account is the number of neighbors, <math display="inline">q</math>, to use in the scheme. In this paper, 8 neighbor points were taken into account following the stencil shown in  [[#img-3|Figure 3]].
 <div id='img-3'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+{| style="text-align: center; border: 1px solid; margin: 1em auto; width: auto;max-width: auto;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig03.png|300x300px|Centered stencil used in the scheme.]]
+|[[Image:Tinoco_Guerrero_et_al_2022b_6735_fig03.png|326px|Centered stencil used in the scheme.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 3:''' Centered stencil used in the scheme.
+| colspan="1" | '''Figure 3'''. Centered stencil used in the scheme
 |}
-Once a discretization of the spatial operator is obtained, the problem of finding an approximation to the solution of the diffusion equation
+Once a discretization of the spatial operator is obtained, following the idea presented in Tinoco-Guerrero et al. <span id='citeF-1'></span>[[#cite-1|[1]]], the semi-discretized PDE can be rewritten as a linear ODE system in time,
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 307: / Line 237: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>    \frac{\partial u}{\partial t} = L_0,   </math>
+| style="text-align: center;" | <math>    \frac{\partial }{\partial t}u_{i}(t) = \sum _{j=0}^{8}\Gamma _{j}u_{j}(t), \qquad  i = 1,\cdots ,m,   </math>
 |}
 |}
-can be solved as a linear system of ordinary differential equations in time; which can be solved through a Runge-Kutta method. In the present work it is proposed to use Runge-Kutta <span id='citeF-8'></span>[[#cite-8|[8]]] of order 2, 3, and 4 as follows:
+where <math display="inline">m</math> is the total amount of grid points.
-==Second order Runge-Kutta==
+There exist several ways to approximate the solution of this ODE system, for example, it is possible to use a forward Euler method as
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 320: / Line 250: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>     u^{k} = u^{k-1} + \frac{\Delta t}{2}\left(r_1 + r_2\right),    </math>
+| style="text-align: center;" | <math>    \frac{u_{i}^{k+1} - u_{i}^{k}}{\Delta t} = \sum _{j=0}^{8}\Gamma _{j}u_{j}^{k},   </math>
 |}
 |}
-where
+nevertheless, this kind of implementation has proven to be ''conditionally stable'' and, in some cases, the stability conditions are difficult to accomplish.
+Another way to solve this system is the Method of Lines, where each <math display="inline">u_i(t)</math> is solved for a fixed grid node <math display="inline">p_i = (x_i, y_i)</math>, ''i.e.'' it is solved by “lines”. As we are now dealing with a set of ODEs, it is possible to use a Runge-Kutta method, to solve the ODE on each line in order to obtain stable and accurate results.
+In the present work it is proposed to use Runge-Kutta <span id='citeF-20'></span>[[#cite-20|[20]]] of 2nd, 3rd, and 4th order, as follows.
+===Second-order Runge-Kutta method===
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 331: / Line 267: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>     r_1 = u^{k-1},    </math>
+| style="text-align: center;" | <math>     u_{i}^{k} = u_{i}^{k-1} + \frac{\Delta t}{2}\left(r_1 + r_2\right),    </math>
 |}
 |}
+where
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 340: / Line 278: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>     r_2 = u^{k-1}+ \frac{\Delta t}{2} (r_1).    </math>
+| style="text-align: center;" | <math>r_1    = u_{i}^{k-1},</math>
+|-
+| style="text-align: center;" | <math>     r_2    = u_{i}^{k-1}+ \frac{\Delta t}{2} (r_1).</math>
+|-
+| style="text-align: center;" | <math>     u_i^{k-1}   = \sum _{j=0}^{8}\Gamma _{j}u_{j}^{k-1}    </math>
 |}
 |}
-==Third order Runge-Kutta==
+===Third-order Runge-Kutta method===
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 351: / Line 293: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>     u^{k} = u^{k-1} + \frac{\Delta t}{6}\left(r_1 + 4r_2 + r_3\right),    </math>
+| style="text-align: center;" | <math>     u_{i}^{k} = u_{i}^{k-1} + \frac{\Delta t}{6}\left(r_1 + 4r_2 + r_3\right),    </math>
 |}
 |}
@@ Line 362: / Line 304: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>     r_1 = u^{k-1},    </math>
+| style="text-align: center;" | <math>r_1    = u_{i}^{k-1},</math>
-|}
-|}
-{| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
-|
+| style="text-align: center;" | <math>     r_2   = u_{i}^{k-1} + \frac{\Delta t}{2} (r_1),</math>
-{| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>     r_2 = u^{k-1} + \frac{\Delta t}{2} (r_1),    </math>
+| style="text-align: center;" | <math>     r_3   = u_{i}^{k-1} + \Delta t (2r_2-r_1),</math>
-|}
-|}
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>     r_3 = u^{k-1} + \Delta t (2r_2-r_1),    </math>
+| style="text-align: center;" | <math>     u_i^{k-1}   = \sum _{j=0}^{8}\Gamma _{j}u_{j}^{k-1}    </math>
 |}
 |}
-==Fourth order Runge-Kutta==
+===Fourth-order Runge-Kutta method===
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 391: / Line 321: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>     u^{k} = u^{k-1} + \frac{\Delta t}{6}\left(r_1 + 2r_2 + 2r_3 + r_4\right),    </math>
+| style="text-align: center;" | <math>     u_{i}^{k} = u_{i}^{k-1} + \frac{\Delta t}{6}\left(r_1 + 2r_2 + 2r_3 + r_4\right),    </math>
 |}
 |}
@@ Line 402: / Line 332: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>     r_1 = u^{k-1},    </math>
+| style="text-align: center;" | <math>r_1    = u_{i}^{k-1},</math>
-|}
-|}
-{| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
-|
+| style="text-align: center;" | <math>     r_2    = u_{i}^{k-1} + \frac{\Delta t}{2} (r_1),</math>
-{| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>     r_2 = u^{k-1} + \frac{\Delta t}{2} (r_1),    </math>
+| style="text-align: center;" | <math>     r_3    = u_{i}^{k-1} + \frac{\Delta t}{2} (r_2),</math>
-|}
-|}
-{| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
-|
+| style="text-align: center;" | <math>     r_4    = u_{i}^{k-1} + \Delta t (r_3),</math>
-{| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>     r_3 = u^{k-1} + \frac{\Delta t}{2} (r_2),    </math>
+| style="text-align: center;" | <math>     u_i^{k-1}   = \sum _{j=0}^{8}\Gamma _{j}u_{j}^{k-1}    </math>
-|}
-|}
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>     r_4 = u^{k-1} + \Delta t (r_3),    </math>
 |}
 |}
@@ Line 435: / Line 346: @@
 in all the cases, <math display="inline">k</math> represents the time level.
-==3 Numerical tests==
+==3. Numerical tests==
-To show the performance of the proposed scheme, the problem of obtaining a numerical solution to diffusion equation on three different regions was proposed. The three regions, denoted as CAB, MID and UP, are non-rectangular planar domains and are shown on Figures [[#img-4a|4a]] to [[#img-5|5]].
+To show the performance of the proposed scheme, the problem of obtaining a numerical solution to the diffusion equation in four different regions was proposed. The first region, denoted as A, corresponds to the unitary square, for comparison reasons, and the other three regions, denoted as B, C, and D, are non-rectangular planar domains. All the regions were scaled to fit on <math display="inline">[0,1]\times [0,1]</math>, and meshed with <math display="inline">21\times{21}</math> nodes, following a variational procedure implemented in UNAMalla <span id='citeF-21'></span>[[#cite-21|[21]]], then they were subdivided to obtain meshes with <math display="inline">41\times{41}</math> nodes. The meshes for with <math display="inline">21\times{21}</math> nodes for each region can be seen in [[#img-4|Figure 4]].
-<div id='img-4a'></div>
-<div id='img-4b'></div>
 <div id='img-4'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig04.png|300x300px|CAB region.]]
+|style="padding:10px;"| [[Image:Tinoco_Guerrero_et_al_2022b_8782_fig04a.png|200x200px]]
-|[[Image:Draft_Tinoco Guerrero_234551248-fig05.png|300x300px|MID region.]]
+|style="padding:10px;"|[[Image:Tinoco_Guerrero_et_al_2022b_4930_fig04b.png|200x200px]]
+|-
+|style="padding:10px;"|[[Image:Tinoco_Guerrero_et_al_2022b_1349_fig04c.png|202x202px]]
+|style="padding:10px;"|[[Image:Tinoco_Guerrero_et_al_2022b_4219_fig04d.png|200x200px|Meshes of the test regions with 21×21 nodes.]]
 |- style="text-align: center; font-size: 75%;"
-| (a) CAB region.
+| colspan="2" style="padding:10px;"| '''Figure 4.''' Meshes of the test regions with <math>21\times{21}</math> nodes
-| (b) MID region.
-|- style="text-align: center; font-size: 75%;"
-| colspan="2" | '''Figure 4'''
 |}
-<div id='img-5'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+For all the regions, two different tests were performed:
+'''Test 1.'''  Following the same idea as in Tinoco-Ruiz <span id='citeF-22'></span>[[#cite-22|[22]]], where a well-known diffusion problem is presented, the initial and boundary conditions were taken from the closed-form solution:
+{| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig06.png|300x300px|UP region.]]
+|
-|- style="text-align: center; font-size: 75%;"
+{| style="text-align: left; margin:auto;width: 100%;"
-| colspan="1" | '''Figure 5:''' UP region.
+|-
+| style="text-align: center;" | <math>     g(x, y, t) = e^{(-2\pi ^2\nu t)}\cos (\pi x)\cos (\pi y).    </math>
+|}
 |}
-All the regions were scaled to fit on <math display="inline">[0,1]\times [0,1]</math>, and meshed with <math display="inline">21\times{21}</math> nodes, following a variational procedure implemented in UNAMalla <span id='citeF-9'></span>[[#cite-9|[9]]], then they were subdivided to obtain meshes with <math display="inline">41\times{41}</math> nodes.
+'''Test 2.'''  For comparison purposes, one of the tests presented in Sánchez et al. <span id='citeF-23'></span>[[#cite-23|[23]]], was also taken into account; in this case, the initial and boundary conditions were taken from the closed-form solution:
-For all the regions, the initial and boundary conditions were taken from the closed form solution:
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 469: / Line 383: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>    g(x, y, t) = e^{(-2\pi ^2\nu t)}\cos (\pi x)\cos (\pi y)   </math>
+| style="text-align: center;" | <math>     g(x, y, t) = e^{(-2\pi ^2\nu t)}\sin (\pi x + \pi y).    </math>
 |}
 |}
-The time interval <math display="inline">[0,1]</math> was subdivided with different discretizations, the number of time steps was chosen to satisfy the classical Courant condition <span id='citeF-10'></span>[[#cite-10|[10]]]:
+The time interval <math display="inline">[0,1]</math> was subdivided with different discretizations, the number of time steps was chosen to satisfy the classical Courant condition <span id='citeF-24'></span>[[#cite-24|[24]]]:
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 486: / Line 400: @@
 where <math display="inline">\Delta x</math> and <math display="inline">\Delta y</math> are the spatial steps on each grid, taken from the meshes with <math display="inline">41\times{41}</math> nodes, and <math display="inline">\nu = 0.2</math>.
-The norm of the quadratic error, at a <math display="inline">k-</math>th time level, can be computed as
+The norm of the quadratic error, at a <math display="inline">k-</math>the time level can be computed as
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 499: / Line 413: @@
 where <math display="inline">U_{i,j}^k</math> and <math display="inline">u_{i,j}^k</math> are the approximated and theoretical solutions, respectively, at the <math display="inline">i,j</math>-th grid node, and <math display="inline">A_{i,j}</math> is the area of the polygon defined by the points <math display="inline">P_{i+1,j}</math>, <math display="inline">P_{i,j+1}</math>, <math display="inline">P_{i-1,j}</math> and <math display="inline">P_{i,j-1}</math>.
-The set of Figures [[#img-6|6]] to [[#img-9|9]] presents a comparison of the numerical results obtained, for some time levels, with the proposed scheme using the fourth order Runge-Kutta approximation and the exact solution for CAB region. The approximation is presenten on the left, and the exact solution on the right.
+[[#img-5|Figure 5]] presents a comparison of the numerical results obtained with the proposed scheme using the second-order Runge-Kutta approximation and the exact solutions for the test region A at time <math display="inline">t = 0.66</math>s. The approximation is presented on the left, and the exact solution on the right.
-<div id='img-6'></div>
+<div id='img-5'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig07.png|500x500px|Comparison between numerical results (left) and exact solution (right) for CAB region at t = 0s.]]
+| [[Image:Review_300670986669-fig13.png|716x716px]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 6:''' Comparison between numerical results (left) and exact solution (right) for CAB region at <math>t = 0</math>s.
-|}
-<div id='img-7'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig08.png|500x500px|Comparison between numerical results (left) and exact solution (right) for CAB region at t = 0.33s.]]
+| [[Image:Review_300670986669-fig14.png|716x716px|Comparison between numerical results (left) and exact solution (right) for test region A. Test 1 (Top) and Test 2 (Bottom).]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 7:''' Comparison between numerical results (left) and exact solution (right) for CAB region at <math>t = 0.33</math>s.
+| style="padding-bottom:10px;"|'''Figure 5'''. Comparison between numerical results (left) and exact solution (right) for test region A. Test 1 (Top) and Test 2 (Bottom)
 |}
-<div id='img-8'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+Similarly, Figures [[#img-6|6]], [[#img-7|7]], and [[#img-8|8]] present the respective comparisons for the regions B, C, and D.
+<div id='img-6'></div>
+{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig09.png|500x500px|Comparison between numerical results (left) and exact solution (right) for CAB region at t = 0.66s.]]
+|[[Image:Review_300670986669-fig15.png|716x716px]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 8:''' Comparison between numerical results (left) and exact solution (right) for CAB region at <math>t = 0.66</math>s.
-|}
-<div id='img-9'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig10.png|500x500px|Comparison between numerical results (left) and exact solution (right) for CAB region at t = 1s.]]
+|[[Image:Review_300670986669-fig16.png|716x716px|Comparison between numerical results (left) and exact solution (right) for test region B. Test 1 (Top) and Test 2 (Bottom).]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 9:''' Comparison between numerical results (left) and exact solution (right) for CAB region at <math>t = 1</math>s.
+| style="padding-bottom:10px;"|'''Figure 6'''. Comparison between numerical results (left) and exact solution (right) for test region B. Test 1 (Top) and Test 2 (Bottom)
 |}
-In a similar way, the sets of Figures [[#img-10|10]] to [[#img-13|13]], and  [[#img-14|14]] to [[#img-17|17]] present the respective comparisons for the MID and the UP regions.
-<div id='img-10'></div>
+<div id='img-7'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig11.png|500x500px|Comparison between numerical results (left) and exact solution (right) for MID region at t = 0s.]]
+|[[Image:Review_300670986669-fig17.png|716x716px]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 10:''' Comparison between numerical results (left) and exact solution (right) for MID region at <math>t = 0</math>s.
-|}
-<div id='img-11'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig12.png|500x500px|Comparison between numerical results (left) and exact solution (right) for MID region at t = 0.33s.]]
+|[[Image:Review_300670986669-fig18.png|716x716px|Comparison between numerical results (left) and exact solution (right) for test region C. Test 1 (Top) and Test 2 (Bottom).]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 11:''' Comparison between numerical results (left) and exact solution (right) for MID region at <math>t = 0.33</math>s.
+| style="padding-bottom:10px;"|'''Figure 7'''. Comparison between numerical results (left) and exact solution (right) for test region C. Test 1 (Top) and Test 2 (Bottom)
 |}
-<div id='img-12'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+<div id='img-8'></div>
+{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig13.png|500x500px|Comparison between numerical results (left) and exact solution (right) for MID region at t = 0.66s.]]
+|[[Image:Review_300670986669-fig19.png|716x716px]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 12:''' Comparison between numerical results (left) and exact solution (right) for MID region at <math>t = 0.66</math>s.
-|}
-<div id='img-13'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig14.png|500x500px|Comparison between numerical results (left) and exact solution (right) for MID region at t = 1s.]]
+|[[Image:Review_300670986669-fig20.png|716x716px|Comparison between numerical results (left) and exact solution (right) for test region D. Test 1 (Top) and Test 2 (Bottom).]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 13:''' Comparison between numerical results (left) and exact solution (right) for MID region at <math>t = 1</math>s.
+| style="padding-bottom:10px;"|'''Figure 8'''. Comparison between numerical results (left) and exact solution (right) for test region D. Test 1 (Top) and Test 2 (Bottom)
 |}
-<div id='img-14'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+Figures [[#img-9|9]] to [[#img-12|12]] present the maximum value of the error <math display="inline">\| e^k\| _2</math> computed for all the regions. For each figure, both tests were performed in the regions meshed with <math display="inline">21\times{21}</math> nodes and <math display="inline">41\times{41}</math> nodes.
+<div id='img-9'></div>
+{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig15.png|500x500px|Comparison between numerical results (left) and exact solution (right) for UP region at t = 0s.]]
+|[[Image:Tinoco_Guerrero_et_al_2022b_3937_fig09a.png|359x359px]]
+|[[Image:Tinoco_Guerrero_et_al_2022b_6678_fig09b.png|359x359px|]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 14:''' Comparison between numerical results (left) and exact solution (right) for UP region at <math>t = 0</math>s.
+| colspan="2" style="padding:10px;"| '''Figure 9'''. Maximum <math>\| e^k\| _2</math> computed for test region A with <math>21\times{21}</math> nodes (Left) and <math>41\times{41}</math> nodes (Right)
 |}
-<div id='img-15'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+<div id='img-10'></div>
+{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig16.png|500x500px|Comparison between numerical results (left) and exact solution (right) for UP region at t = 0.33s.]]
+|[[Image:Tinoco_Guerrero_et_al_2022b_3582_fig10a.png|359px]]
+|[[Image:Tinoco_Guerrero_et_al_2022b_8582_fig10b.png|359px]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 15:''' Comparison between numerical results (left) and exact solution (right) for UP region at <math>t = 0.33</math>s.
+| colspan="2" style="padding:10px;"| '''Figure 10'''. Maximum <math>\| e^k\| _2</math> computed for test region B with <math>21\times{21}</math> nodes (Left) and <math>41\times{41}</math> nodes (Right)
 |}
-<div id='img-16'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+<div id='img-11'></div>
+{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig17.png|500x500px|Comparison between numerical results (left) and exact solution (right) for UP region at t = 0.66s.]]
+|[[Image:Tinoco_Guerrero_et_al_2022b_6566_fig11a.png|359px]]
+|[[Image:Tinoco_Guerrero_et_al_2022b_7935_fig11b.png|359px]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 16:''' Comparison between numerical results (left) and exact solution (right) for UP region at <math>t = 0.66</math>s.
+| colspan="2" style="padding:10px;"| '''Figure 11'''. Maximum <math>\| e^k\| _2</math> computed for test region C with <math>21\times{21}</math> nodes (Left) and <math>41\times{41}</math> nodes (Right)
 |}
-<div id='img-17'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+<div id='img-12'></div>
+{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
 |-
-|[[Image:Draft_Tinoco Guerrero_234551248-fig18.png|500x500px|Comparison between numerical results (left) and exact solution (right) for UP region at t = 1s.]]
+|[[Image:Tinoco_Guerrero_et_al_2022b_8962_fig12a.png|359px]]
+|[[Image:Tinoco_Guerrero_et_al_2022b_2351_fig12b.png|359px]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 17:''' Comparison between numerical results (left) and exact solution (right) for UP region at <math>t = 1</math>s.
+| colspan="2" style="padding:10px;"| '''Figure 12'''. Maximum <math>\| e^k\| _2</math> computed for test region D with <math>21\times{21}</math> nodes (Left) and <math>41\times{41}</math> nodes (Right)
 |}
-Table [[#table-1|1]] shows the results of the computed quadratic errores computed for all the test on CAB region. All the methods are reported with the following names: RK''N'', where ''N'' represents the Runge-Kutta order. Analogously, Tables [[#table-2|2]] and [[#table-3|3]] present the results for MID and UP regions respectively.
+==4. Conclusions==
+The numerical results show that the proposed method of Lines applied to the diffusion equation produces satisfactory numerical solutions. In the tests carried out, no spurious oscillations or instabilities were perceived. In addition, the results show that it is not necessary to “work more” by using higher Runge-Kutta methods; this scheme accomplishes better results when second-order Runge-Kutta is used to solve the Ordinary Differential Equations system.
-{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
+It is possible to appreciate that with the presented discretizations, even with the fourth-order Runge-Kutta method, acceptable numerical results can be obtained. Additionally, the numerical results shown in Figures [[#img-9|9]] to [[#img-12|12]] show that a refinement of the spatial mesh can improve the approximations carried out with the method; as expected in these kinds of methods. It is worth mentioning that, for all the tests, a minimum number of time steps was used, taking into account the classical theory for the diffusion equation, which makes the computational cost of this method low, since it does not require making very small temporary discretizations, as in other cases.
-|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 Computed <math>\| e^k\| _2</math> for CAB.
-|- style="border-top: 2px solid;"
-| colspan='4' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | CAB region
-|- style="border-top: 2px solid;"
-| style="border-left: 2px solid;border-right: 2px solid;" |             Method
-| style="border-left: 2px solid;border-right: 2px solid;" | time
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>21\times{21}</math>nodes
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>41\times{41} nodes</math>
-|- style="border-top: 2px solid;"
-| rowspan='4' style="border-left: 2px solid;border-right: 2px solid;" | RK2
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.25</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.1987E-04</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>2.6751E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.50</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>4.5951E-05</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.0264E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.75</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.7339E-05</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>3.8606E-06</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.00</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>6.5273E-06</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.4495E-06</math>
-|- style="border-top: 2px solid;"
-| rowspan='4' style="border-left: 2px solid;border-right: 2px solid;" | RK3
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.25</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.2515E-04</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>3.172E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.50</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>4.7959E-05</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.2136E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.75</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.8110E-05</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>4.5797E-06</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.00</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>6.8218E-06</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.7241E-06</math>
-|- style="border-top: 2px solid;"
-| rowspan='4' style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | RK4
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.25</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>4.7901E-05</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.2122E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.50</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>6.7934E-06</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.7169E-06</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.75</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>9.5131E-07</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>2.4025E-07</math>
-|- style="border-bottom: 2px solid;"
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.00</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.3353E-07</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>3.3712E-08</math>
-|}
+Furthermore, the proposed scheme can produce stable results for non-standard initial conditions in highly irregular regions. For example, in the following videos
+* https://youtu.be/soUxe5uUY_U
+* https://youtu.be/479YeBjqmSs
-{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
-|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 Computed <math>\| e^k\| _2</math> for MID.
-|- style="border-top: 2px solid;"
-| colspan='4' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | MID region
-|- style="border-top: 2px solid;"
-| style="border-left: 2px solid;border-right: 2px solid;" |             Method
-| style="border-left: 2px solid;border-right: 2px solid;" | time
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>21\times{21}</math>nodes
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>41\times{41} nodes</math>
-|- style="border-top: 2px solid;"
-| rowspan='4' style="border-left: 2px solid;border-right: 2px solid;" | RK2
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.25</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>6.0947E-04</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.3452E-04</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.50</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>2.3585E-04</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>5.1647E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.75</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>8.8078E-05</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.9274E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.00</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>3.2834E-05</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>7.1847E-06</math>
-|- style="border-top: 2px solid;"
-| rowspan='4' style="border-left: 2px solid;border-right: 2px solid;" | RK3
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.25</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>6.0920E-04</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.3433E-04</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.50</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>2.3574E-04</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>5.1570E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.75</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>8.8040E-05</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.9245E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.00</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>3.2819E-05</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>7.1739E-06</math>
-|- style="border-top: 2px solid;"
-| rowspan='4' style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | RK4
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.25</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>2.3575E-04</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>5.1571E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.50</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>3.2817E-05</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>7.1734E-06</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.75</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>4.5586E-06</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>9.9644E-07</math>
-|- style="border-bottom: 2px solid;"
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.00</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>6.3330E-07</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.3843E-07</math>
-|}
+solutions of the diffusion equation are presented on domains that are geometrical approximations of real geographical locations, where the boundary conditions are fixed as <math display="inline">0</math> and the initial condition is stated as:
+{| class="formulaSCP" style="width: 100%; text-align: left;"
-{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
-|+ style="font-size: 75%;" |<span id='table-3'></span>Table. 3 Computed <math>\| e^k\| _2</math> for UP.
-|- style="border-top: 2px solid;"
-| colspan='4' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | UP region
-|- style="border-top: 2px solid;"
-| style="border-left: 2px solid;border-right: 2px solid;" |             Method
-| style="border-left: 2px solid;border-right: 2px solid;" | time
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>21\times{21}</math>nodes
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>41\times{41} nodes</math>
-|- style="border-top: 2px solid;"
-| rowspan='4' style="border-left: 2px solid;border-right: 2px solid;" | RK2
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.25</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>7.7062E-04</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>2.0092E-04</math>
 |-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.50</math>s
+|
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>3.8002E-04</math>
+{| style="text-align: left; margin:auto;width: 100%;"
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>9.7678E-05</math>
 |-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.75</math>s
+| style="text-align: center;" | <math>    g(x,y) =  \begin{cases}       1 & x,y \hbox{ in } [0.35, 0.65]\\       0 & \hbox{ All the other cases.}      \end{cases}   </math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.5808E-04</math>
+|}
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>4.0295E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.00</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>6.1880E-05</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.5696E-05</math>
-|- style="border-top: 2px solid;"
-| rowspan='4' style="border-left: 2px solid;border-right: 2px solid;" | RK3
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.25</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>7.7432E-04</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>2.0465E-04</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.50</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>3.8204E-04</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>9.9664E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.75</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.5896E-04</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>4.1144E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.00</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>6.2233E-05</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.6032E-05</math>
-|- style="border-top: 2px solid;"
-| rowspan='4' style="border-bottom: 2px solid;border-left: 2px solid;border-right: 2px solid;" | RK4
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.25</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>3.8188E-04</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>9.9622E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.50</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>6.2129E-05</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.6005E-05</math>
-|-
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.75</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>8.9096E-06</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>2.2853E-06</math>
-|- style="border-bottom: 2px solid;"
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.00</math>s
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.2497E-06</math>
-| style="border-left: 2px solid;border-right: 2px solid;" | <math>3.2014E-07</math>
 |}
-==4 Conclusions==
+It is possible to see that even for these conditions in these regions, the scheme produces stable results that show the expected behavior.
-It can be observed in the results presented in this work that the Line Method proposed for the solution of the diffusion equation produces satisfactory numerical solutions of the equation. In the tests carried out, no spurious oscillations or instabilities were perceived; Furthermore, it is shown that the approximations are more precise when using a fourth-order Runge-Kutta method.
+An important remark is that, even when in this work structured convex grids were used, the development of the method doesn't take into account a particular data structure, i.e. this method can be used not only on structured meshes but also as a meshless method, that would be developed as future work.
-It is possible to appreciate that with the presented discretizations, even with the second-order Runge-Kutta method, acceptable numerical results can be obtained. Additionally, the numerical results in Tables [[#table-1|1]] to [[#table-3|3]] show that a refinement of the spatial mesh can improve the approximations carried out with the method; as expected in this kind of methods. It is worth mentioning that, for all the tests, a minimum number of time steps was used, taking into account the classical theory for the diffusion equation, which makes the computational cost of this method low, since it does not require making very small temporary discretizations, as in other cases.
+==Acknowledgments==
-==5 Acknowledgments==
+We want to thank AULA CIMNE-Morelia and CIC-UMSNH for the financial support for this work. Thanks to Universidad Vasco de Quiroga for lending the facilities to perform the final part of the work done for the paper.
-We want to thank AULA CIMNE-Morelia and CIC-UMSNH for the financial support for this work.
 Special thanks to the ``Laboratorio de Matemáticas Aplicadas'' of the ``Facultad de Ciencias Físico-Matemáticas'' at ``Universidad Michoacana de San Nicolás de Hidalgo'' for lending their computational infrastructure.
-===BIBLIOGRAPHY===
+The authors would also like to thank the reviewers for their valuable feedback and suggestions to improve the present work.
+==References==
+<div class="auto" style="text-align: left;width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">
 <div id="cite-1"></div>
-'''[[#citeF-1|[1]]]''' Gerardo Tinoco-Guerrero and Francisco Javier Domínguez-Mota and José Gerardo Tinoco-Ruiz and Juan Salvador Lucas-Martínez. (2018) "Aproximación de la Ecuación de Advección en Regiones Irregulares Utilizando un Método de Líneas y Diferencias Finitas Generalizadas", Volume 2. Revista Mexicana de Métodos Numéricos 2
+[[#citeF-1|[1]]]  Tinoco-Guerrero G.,  Domínguez-Mota F.J., Tinoco-Ruiz J.G., Lucas-Martínez J.S.  Aproximación de la ecuación de advección en regiones irregulares utilizando un método de líneas y diferencias finitas generalizadas. Revista Mexicana de Métodos Numéricos, Volume 2, 3, 2028.
 <div id="cite-2"></div>
-'''[[#citeF-2|[2]]]''' Victorita Dolean and Pierre Jolivet and Frédéric Nataf. (2015) "An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation". Society for Industrial and Applied Mathematics
+[[#citeF-2|[2]]] Manshoor B.,  Salleh H., Khalid A., Sayed Abdelaal M.A. Method of lines and Runge-Kutta method in solving partial differential equation for eeat equation. Journal of Complex Flow, 3(1):21-25, 2021.
 <div id="cite-3"></div>
-'''[[#citeF-3|[3]]]''' William Schiesser. (1991) "The Numerical Method of Lines". San Diego : Academic Press
+[[#citeF-3|[3]]] Ureña Prieto F., Benito Muñoz J.J.,  Gavete Corvinos L. Application of the generalized finite difference method to solve the advection-diffusion equation. Journal of Computational and Applied Mathematics, 25:1849-1855, 2011.
 <div id="cite-4"></div>
-'''[[#citeF-4|[4]]]''' Gerardo Tinoco-Guerrero and Francisco Javier Domínguez-Mota and José Gerardo Tinoco-Ruiz. (2015) "Numerical Solution of Advection Equations Using Structured Grids on Plane Irregular Regions With a Finite Differences Scheme", Volume 1. Boletón de la Sociedad Mexicana de Computación Científica y Sus Aplicaciones 3 29 - 32
+[[#citeF-4|[4]]] Li P.-W.,  Fan C.-M., Grabski J.K. A meshless generalized finite difference method for solving shallow water equations with the flux limiter technique. Engineering Analysis with Boundary Elements, 131:159&#8211;173, 2021.
 <div id="cite-5"></div>
-'''[[#citeF-5|[5]]]''' John C. Strikwerda. (2004) "Finite Difference Schemes and Partial Differential Equations". Society for Industrial and Applied Mathematics
+[[#citeF-5|[5]]] Zhang T.,  Lin Z.-H.,  Huang G.-Y.,  Fan C.-M., Li P.-W. Solving Boussinesq equations with a meshless finite difference method. Ocean Engineering, 198:106957, 2020.
 <div id="cite-6"></div>
-'''[[#citeF-6|[6]]]''' J. W. Thomas. (1995) "Numerical Partial Differential Equations: Finite Difference Methods", Volume 2. Springer
+[[#citeF-6|[6]]] Fu Z.-J.,  Tang Z.-C.,  Zhao H.-T., Li P.-W.,  Rabczuk T. Numerical solutions of the coupled unsteady nonlinear convection-diffusion equations based on generalized finite difference method. The European Physical Journal Plus, 134(6):1&#8211;20, 2019.
 <div id="cite-7"></div>
-'''[[#citeF-7|[7]]]''' Gerardo Tinoco-Guerrero and Francisco Javier Domínguez-Mota and José Gerardo Tinoco-Ruiz. (2020) "A Study of the Stability for a Generalized Finite-Difference Scheme Applied to the Advection-Diffusion Equation", Volume 176. Mathematics and Computer in Simulations 301 - 311
+[[#citeF-7|[7]]] Zhang T.,  Ren Y.-F.,  Yang Z.-Q.,  Fan C.-M., Li P.-W. Application of generalized finite difference method to propagation of nonlinear water waves in numerical wave flume. Ocean Engineering, 123:278&#8211;290, 2016.
 <div id="cite-8"></div>
-'''[[#citeF-8|[8]]]''' J. C. Butcher. (1969) "Conference on the Numerical Solution of Differential Equations", Volume 109. Springer 133 - 139
+[[#citeF-8|[8]]] Li P.-W., Fan C.-M. Generalized finite difference method for two-dimensional shallow water equations. Engineering Analysis with Boundary Elements, 80:58&#8211;71, 2017.
 <div id="cite-9"></div>
-'''[[#citeF-9|[9]]]''' UNAMalla. (2011) "An Automatic Package for Numerical Grid Generation"
+[[#citeF-9|[9]]] Fu Z.-J., Zhang J., Li P.-W.,  Zheng J.-H. A semi-Lagrangian meshless framework for numerical solutions of two-dimensional sloshing phenomenon. Engineering Analysis with Boundary Elements, 112:58&#8211;67, 2020.
 <div id="cite-10"></div>
-'''[[#citeF-10|[10]]]''' Michael B. Abbott and Andrew D. McCowan and Ian R. Warren. (1984) "Accuracy of Short-Wave Numerical Models", Volume 110. Journal of Hydraulic Engineering 10 1287 - 1301
+[[#citeF-10|[10]]] Wang Q., Kim P., Qu W. A hybrid localized meshless method for the solution of transient groundwater flow in two dimensions. Mathematics, 10(3):1-14, 2022.
+<div id="cite-11"></div>
+[[#citeF-11|[11]]] Domínguez-Mota F.J.,  Mendoza-Armenta S.,  Tinoco-Ruiz J.G. Finite difference schemes satisfying an optimality condition. IMACS Series in Computational and Applied Mathematics: MASCOT, 2011.
+<div id="cite-12"></div>
+[[#citeF-12|[12]]] Tinoco-Guerrero G., Domínguez-Mota F.J.,  Tinoco-Ruiz J.G. Numerical solution of advection equations using structured grids on plane irregular regions with a finite differences scheme. Boletín de la Sociedad Mexicana de Computación Científica y Sus Aplicaciones, 1(3):29-32, 2015.
+<div id="cite-13"></div>
+[[#citeF-13|[13]]] Tinoco-Ruiz J.G.,  Domínguez-Mota F.J.,  Tinoco-Guerrero G., Fernández-Valdés P.M., Mendoza-Armenta S. Numerical solution of differential equations in irregular plane regions using quality structured convex grids.  International Journal of Modeling, Simulation and Scientific Computing, 4(2):1340005, 2013.
+<div id="cite-14"></div>
+[[#citeF-14|[14]]] Domínguez-Mota F.J., Lucas-Martínez J.S., Tinoco-Guerrero G. A generalized finite-differences scheme used in modeling of a direct and an inverse problem of advection-diffusion. International Journal of Applied Mathematics, 33(4):599-608, 2020.
+<div id="cite-15"></div>
+[[#citeF-15|[15]]] Chávez-Negrete C.,  Domínguez-Mota F.J.,  Santana-Quinteros D. Numerical solution of Richards' equation of water flow by generalized finite differences. Computers and Geotechnics, 101:168&#8211;175, 2018.
+<div id="cite-16"></div>
+[[#citeF-16|[16]]] Schiesser W. The numerical method of lines. San Diego: Academic Press, 1991.
+<div id="cite-17"></div>
+[[#citeF-17|[17]]] Strikwerda J.C. Finite difference schemes and partial differential equations. Society for Industrial and Applied Mathematics, pp. 427, 2004.
+<div id="cite-18"></div>
+[[#citeF-18|[18]]] Thomas J.W. Numerical partial differential equations: finite difference methods. Texts in Applied Mathematics (TAM), Volume 22, Springer, 1995.
+<div id="cite-19"></div>
+[[#citeF-19|[19]]] Tinoco-Guerrero G., Domínguez-Mota F.J.,  Tinoco-Ruiz J.G. A study of the stability for a generalized finite-difference scheme applied to the advection-diffusion equation. Mathematics and Computer in Simulations, 176:301-311, 2020.
+<div id="cite-20"></div>
+[[#citeF-20|[20]]]  Butcher J.C. Conference on the numerical solution of differential equations.  Lecture Notes in Mathematics (LNM), Springer, 109:133-139, 1969.
+<div id="cite-21"></div>
+[[#citeF-21|[21]]] UNAMalla. An automatic package for numerical grid generation, 2011.
+<div id="cite-22"></div>
+[[#citeF-22|[22]]] Tinoco-Ruiz J.G.,  Domínguez-Mota F.J.,  Tinoco-Guerrero G., Lucas-Martínez J.S. Study of the stability of a generalized finite difference scheme applied to the diffusion equation in irregular 2D space regions using convex grids. IMACS Series in Computational and Applied Mathematics: MASCOT2018, 2018.
+<div id="cite-23"></div>
+[[#citeF-23|[23]]]  Sánchez L.M.,  Ureña F., Benito J.J.,  Gavete L. Resolución de la ecuación de difusión en 2-D y 3-D utilizando diferencias finitas generalizadas. Consistencia y Estabilidad. XXI Congreso de Ecuaciones Diferenciales y Aplicaciones,  Ciudad Real, pp. 1-8, 2009.
+<div id="cite-24"></div>
+[[#citeF-24|[24]]]  Abbott M.B.,  McCowan A.D.,  Warren I.R. Accuracy of short-wave numerical models. Journal of Hydraulic Engineering, 110(10):1287-1301, 1984.