## 4 Numerical Experiments  
In this section, we reproduce nontrivial analytical solutions in one space dimension. The Dirichlet boundary conditions at the two end points of the interval and the initial \(u^{0}(x,\ t)\) function were obtained by simply substituting the appropriate \(x\) and \(t\) values into the analytical solution. The MATLAB R2019b software was used for all numerical calculations. We calculated the eigenvalues of the \(M\) matrix in Equation (16) for this problem. Using these, the stiffness ratio and the largest time step size where the standard FTCS method is stable can be determined as usual [33]. To measure the accuracy, the widely used \(L_{\infty}\) error is calculated, which compares the exact value \(u^{\text{analytic}}_{i}\) and the result \(u^{\text{num}}_{i}\) obtained by the actual numerical method at the final time \(t^{\text{fin}}\):  
\[Error=\max_{1\leq i\leq N}\Big{|}u^{\text{analytic}}_{i}(t^{\text{fin}})-u^{ \text{num}}_{i}(t^{\text{fin}})\Big{|}. \tag{19}\]  
Experiment 1. We are going to reproduce the following reference solution, which is a quite recent [48] analytical solution:  
\[u^{\text{analytic}}(x,t)=\frac{x}{t^{7/2}}\bigg{(}1-\frac{x^{2}}{3\alpha t}+ \frac{x^{4}}{60\alpha^{2}t^{2}}\bigg{)}e^{-\frac{x^{2}}{64\alpha t}}. \tag{20}\]  
We considered the domain with the boundaries \(x_{0}=-5\), \(x_{N}=5\), \(t^{0}=0.3\), \(t^{\text{fin}}=0.56\), while \(\alpha=1\). We set \(\Delta x=0.02\), which implies that the stiffness ratio is \(1.01\cdot 10^{5}\), while \(\Delta t^{\text{FTCS}}_{\text{MAX}}=3.9\cdot 10^{-4}\). We examined the error defined in (19) as a function of the time step size \(\Delta t\). First, the errors were calculated for a very large \(\Delta t\), then with decreased time step sizes until small error values were obtained. In Figures 2 and 3, the errors as a function of the time step size are shown in log-log diagrams. Since we used a fixed-space step size and decreased only the time step size, the errors cannot go to zero, but only to the residual error due to space discretization, which is given by (15). This can be seen in the bottom left of the figures. One can notice that although the CLL method is third-order and it is much more accurate than the LNe method, it converges more slowly than the hopscotch-type methods. In accordance with some experiments in our previous papers, the LH method was the most accurate. Figures 4-6 show how the numerical solutions approximate the analytical solution as the time step size is decreased. In the case of the one-stage CNe method and the three-stage CLL methods, the curves are smooth without unphysical oscillations, while the ASH method produces significant oscillations for larger time step sizes. On the other hand, the CNe method converges very slowly and the CLL converges at a medium rate, while the ASH converges much faster. Some investigation of the behaviour of the errors are presented in our previous paper [39], as well for fixed time step size and space-dependent diffusion coefficient. Figure 3: The errors as a function of the time step size for the quickly converging methods, as well as the CLL method in the case of Experiment 1. The error curves reach the residual error \(\tau_{0}\) on the left side of the figure. Figure 2: Errors as a function of the time step size for the slowly converging methods. Figure 4: The numerically calculated temperatures \(u\) for the three-stage CNe method as a function of the space variable \(x\) for different time step sizes in Experiment 1. The red and the black solid lines are for the initial and the analytically calculated temperatures, respectively. Figure 5: The numerically calculated temperatures \(u\) for the three-stage CLL method as a function of the space variable \(x\) for different time step sizes in Experiment 1. The red and the black solid lines are for the initial and the analytically calculated temperatures, respectively. Experiment 2.