## 2 The Analyzed Numerical Methods  
We considered the 1D interval \(x\in\left[x_{0}\,,\,\,x_{N}\right]\) with length \(L=x_{N}-x_{0}\) and constructed an equidistant spatial grid with the node coordinates \(x_{0}\,,\,x_{1}\,,\,\ldots\,,\,x_{N}\), so \(x_{j}=x_{j-1}+\Delta x\), \(j=1,\ldots,N\), \(\Delta x=L/N\). The time variable is also discretized, so if \(t\in\left[t^{0},\,t^{\text{fin}}\right]\), then \(t^{j}=t^{0}+j\Delta t\,,\,\,j=1,\ldots,\,T,T\Delta t=t^{\text{fin}}-t^{0}\). The mesh-ratio is defined as usual, \(r=\frac{a\,\Delta t}{\Delta x^{2}}\). In this work, only the simplest case of Equation (1) (one space-dimensional equidistant mesh and constant \(\alpha\) coefficient) was examined, but we give the original publications where more details can be found about the application of the methods for general meshes as well. Let us start with a brief presentation of the defining formulas of the methods. The UPFD method is proposed by [PERSON] and [PERSON] [41] for the linear diffusion-advection-reaction PDE. It is a non-standard combination of the explicit and implicit Euler discretizations, where only the actual node is considered implicitly and the neighbors explicitly. In the case of Equation (1) it reads: \[\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}=\frac{\alpha}{\Delta x^{2}}\Big{(}u_{i-1 }^{n}-2u_{i}^{n+1}+u_{i+1}^{n}\Big{)}.\] (2) The fully explicit form is: \[u_{i}^{n+1}=\frac{u_{i}^{n}+r\big{(}u_{i-1}^{n}+u_{i+1}^{n}\big{)}}{1+2r},\] (3)
2. The first scheme we invented is the so-called constant neighbor (CNe) algorithm [32; 42], which has the following formula: \[u_{i}^{n+1}=u_{i}^{n}\cdot e^{-2r}+\frac{u_{i-1}^{n}+u_{i+1}^{n}}{2}\Big{(}1-e ^{-2r}\Big{)},\] (4)
3. The CpC method [31] is the organization of the CNe scheme into a two-stage algorithm. The first stage is a fractional time step with length \(p\Delta t\) with the CNe formula: \[u_{i}^{\text{Cp}}=u_{i}^{n}e^{-2 pr}+\frac{u_{i-1}^{n}+u_{i+1}^{n}}{2}\Big{(}1- e^{-2 pr}\Big{)}.\] These new and the old values are combined linearly to obtain the predictor values: \[u_{i}^{\text{pred}}=\Big{(}1-\frac{1}{2p}\Big{)}u_{i}^{n}+\frac{1}{2p}u_{i}^{ \text{Cp}}.\] These are used at the second stage, which is a full-time step size corrector step; thus, the final values are: \[u_{i}^{n+1}=u_{i}^{n}\cdot e^{-2r}+\frac{u_{i-1}^{\text{pred}}+u_{i+1}^{\text{ pred}}}{2}\Big{(}1-e^{-2r}\Big{)},\] (5)
4. CpCC: We now make an attempt to improve the accuracy of the CpC method by iterating one more time using the values obtained by the CpC method as predictors. So, after executing Formula (5) in the CpC method, we set \(u_{i}^{\text{pred}}=u_{i}^{n+1}\) and then calculate Formula (5) again. 5. The two-stage linear-neighbor (LNe or Lne2) method [32] starts with taking a full-size time step with the CNe method to obtain the predictor values \(u_{i}^{\text{pred}}\), which are valid at the end of the current time step. Using these values, a quantity can be introduced: \[s_{i}=u_{i-1}^{\text{pred}}+u_{i+1}^{\text{pred}}-u_{i-1}^{n}-u_{i+1}^{n},\] which is proportional to the aggregated or effective slopes of the neighbors of node \(i\). Now the corrector values of the two-stage LNe method are: \[u_{i}^{n+1}=u_{i}^{n}e^{-2r}+\frac{u_{i-1}^{n}+u_{i+1}^{n}}{2}\Big{(}1-e^{-2r }\Big{)}+\frac{s_{i}}{2}\bigg{(}1-\frac{1-e^{-2r}}{2r}\bigg{)}.\] (6)
6. The LNe3 method is a three-stage algorithm [32] whose first two stages are the same as the LNe schemes. At the end of its second stage, based on the corrector values in Equation (6), one can set \(u_{i}^{\text{pred}}=u_{i}^{n+1}\), recalculate the slopes \(s_{i}\) and then repeat (6) to obtain new corrector results. This procedure gives a three-stage scheme altogether. This algorithm is still second-order, but more accurate than LNe2. The LNe4 method is a four-stage algorithm, which is obtained by repeating the procedure explained in the previous point after the calculations of the LNe3 method.