## 5 Discussion and Summary  
We have studied 20 numerical methods for the non-steady-state linear diffusion equation. All of the algorithms are unconditionally stable explicit methods; 16 of them have been recently constructed by our research group. We have calculated the truncation errors for each algorithm, which was unknown until now in several cases. Moreover, we have considered another method of investigation, in which the numerical solution is compared to the concrete analytical solution of the ODE system, which is obtained by the spatial discretization of the PDE. We think that this latter method produces errors which explain the observed numerical order of the algorithms more clearly if the spatial discretization is kept fixed and only the time step size is decreased, especially when the stiffness ratio is high. We have reproduced a recent analytical solution by the studied 20 numerical methods in six numerical experiments. We have made attempts to explain the observed numerical errors based on the analytically calculated truncation errors. We emphasize again that although these algorithms are only conditionally consistent, they give accurate results for time step sizes orders of magnitude larger than the standard explicit methods, and faster (for large systems) than the implicit methods. Conceptualization, supervision and resources, E.K. ; methodology, E.K., A.N. and J.M. ; analytical investigation, J.M. and A.N. ; software, numerical investigation and visualization, A.N. ; writing--original draft preparation, E.K. and J.M. ; writing--review and editing, A.N. All authors have read and agreed to the published version of the manuscript. This research received no external funding. Here, we present those errors which we were able to calculate. 1. UPFD  
\[\tau_{\text{UPDF}}=\tau_{0}+2\alpha^{2}\frac{u_{(2x)}}{\Delta x^{2}}\Delta t +\frac{\alpha^{2}}{2}u_{(4x)}\Delta t+\alpha^{3}u_{(4x)}\frac{\Delta t^{2}}{ \Delta x^{2}}+\frac{\alpha^{3}}{6}u_{(6x)}\Delta t^{2}\]  
\[\varepsilon_{\text{UPDF}}=\varepsilon_{0}+2\alpha^{2}\frac{u_{(2x)}}{\Delta x ^{2}}\Delta t+\frac{2}{3}\alpha^{2}u_{(4x)}\Delta t+\frac{\alpha^{2}}{180}u_{ (6x)}\Delta x^{2}\Delta t-4\alpha^{3}u_{(2x)}\frac{\Delta t^{2}}{\Delta x^{4} }-\frac{\alpha^{3}}{3}u_{(4x)}\frac{\Delta t^{2}}{\Delta x^{2}}+\frac{7}{45} \alpha^{3}u_{(6x)}\Delta t^{2}\]  
\[\delta_{\text{UPFD}}=-\frac{\beta^{2}}{2}\big{(}u_{i\pm 1}^{n}-2u_{i}^{n} \big{)}\Delta t+\frac{5\beta^{3}}{2}\big{(}u_{i\pm 1}^{n}-2u_{i}^{n}\big{)} \Delta t^{2}+O\Big{(}\Delta t^{3}\Big{)}.\]  
2. CNe  
\[\tau_{\text{CNe}}=\varepsilon_{\text{CNe}}=\varepsilon_{0}+\alpha^{2}\frac{u _{(2x)}}{\Delta x^{2}}\Delta t+\frac{7}{12}\alpha^{2}u_{(4x)}\Delta t+\frac{ \alpha^{2}}{360}u_{(6x)}\Delta x^{2}\Delta t-\frac{2}{3}\alpha^{3}u_{(2x)} \frac{\Delta t^{2}}{\Delta x^{4}}-\frac{\alpha^{3}}{18}u_{(4x)}\frac{\Delta t ^{2}}{\Delta x^{2}}+\frac{89}{540}\alpha^{3}u_{(6x)}\Delta t^{2}\]  
\[\delta_{\text{CNe}}=\frac{\beta^{2}}{2}\big{(}u_{i\pm 1}^{n}-2u_{i}^{n}\big{)} \Delta t-\frac{5\beta^{3}}{2}\big{(}u_{i\pm 1}^{n}-2u_{i}^{n}\big{)} \Delta t^{2}+O\Big{(}\Delta t^{3}\Big{)}\]  
3. CpC. If the parameter \(p\) is kept, we have:  
\[\varepsilon_{\text{CpC}}=\varepsilon_{0}+\frac{\alpha^{2}}{12}u_{(6x)} \Delta x^{2}\Delta t+\frac{\alpha^{3}\Big{(}360(1+3p)u_{(2x)}+30(19+21p) \Delta x^{2}u_{(4x)}+(271+93p)\Delta x^{4}u_{(6x)}\Big{)}}{1080\Delta x^{4}} \Delta t^{2}.\]One can see that the error decreases with decreasing \(p\). However, the stability analysis in our original paper showed that \(p\) must be at least \(\frac{1}{2}\) for unconditional stability, thus this is the optimal choice for \(p\).