## 1 Introduction  
Fickian diffusion or Fourier-type heat conduction is one of the most fundamental mass- or energy-transport processes. In the simplest case, it can be described by a single linear partial differential equation (PDE) of space and time:  
\[\frac{\partial u(x,t)}{\partial t}=\alpha\frac{\partial^{2}u(x,t)}{\partial x^ {2}} \tag{1}\]  
where \(x,\ t\in\mathbb{R}\), \(u=u(x,t)\) is the concentration (temperature in case of heat conduction) as a function of space and time, and \(\alpha\) is the constant diffusion coefficient. The generalizations of this equation, such as the advection (or drift)-diffusion equation and different kinds of advection-diffusion-reaction equations [1], can model the transport and propagation of various particles in physical, chemical and biological systems [2]. Moreover, very similar, albeit more complicated, PDEs are used to model fluid flow through porous media [3], such as moisture [4], ground water or crude oil in underground reservoirs [5]. It is well-known that countless numerical methods have been proposed to solve Equation (1) and its generalizations. Most of them belong to the broad family of finite difference schemes (FDM) [6, 7], which often means a kind of method of lines [1], where the spatial variables are discretized first to obtain an ordinary differential equation (ODE) system and then an ODE solver is employed. Nevertheless, finite element methods (FEM) are also used [8], especially if the geometry is complex. Most of these methods can be considered either as an explicit or an implicit method, and sometimes they have a mixed character [9]. If the physical properties, i.e., the coefficients of the equation, are highly non-uniform in space, the eigenvalues of the system matrix likely have a range of several orders of magnitude, which implies that one has a highly stiff problem. In this case, the so-called Courant-Friedrichs-Lewy (CFL) limit (which refers to conventional explicit methods, e.g., the Runge-Kutta or Adams-Bashforth schemes) is very small. This means that these algorithms are unstable unless the time step size is smaller than this low threshold, and therefore the simulation can be unacceptably slow. Since implicit methods have much better stability properties, they are frequently used for these equations [10; 11; 12; 13; 14; 15; 16]. Nevertheless, they involve the solution of a system of algebraic equations at each time step whose parallelization is much less straightforward than that of explicit methods. Thus, the calculations can be time-consuming for very-large-sized and non-tridiagonal matrices, which is the case when the number of space dimensions is larger than one. In these cases, even the simplest explicit Euler time integration can be much faster than the implicit one [10; 11; 17]. Regrettably, the formerly rapid increase of the clock frequencies of processors disappeared about two decades ago, which reinforced the progression towards increasing parallelism in high-performance computing [18; 19]. These are the main reasons why there is a growing interest in explicit and unconditionally stable algorithms [20; 21; 22; 23; 24; 25; 26]. In the last few years, our research group has developed a couple of new algorithms of this kind. In our original papers [14; 27; 28; 29; 30; 31; 32; 33], we examined our new methods analytically and proved that they are indeed unconditionally stable for the linear heat equation. However, when we examined the consistency, we always found that they were only conditionally consistent, in accordance with the literature [11; 12; 25; 34] (p. 120). The goal of this paper is to collect a large number of these methods, give their truncation error, and draw some consequences about the accuracy of the methods. Note that in textbooks and monographs, typically the standard methods, such as the explicit and implicit Euler, as well as the Crank-Nicolson method, are analyzed [35; 36].