## 1 Introduction  
Diffusion of particles and Fourier-type heat conduction are omnipresent mass or energy transport processes. In the simplest linear case, they are described by the following partial differential equation (PDE):  
\[\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}\) are the independent variables, \(u=u(x,t)\) is the unknown concentration of particles or the temperature in the case of heat transfer, and \(\alpha\) is the coefficient of (thermal) diffusivity. The thermal diffusivity of a material can be given as \(\alpha=k/(c\rho)\), where \(c=c(\overset{\rightarrow}{r},t)\), \(k=k(\overset{\rightarrow}{r},t)\), and \(\rho=\rho(\overset{\rightarrow}{r},t)\) are the specific heat, the heat conductivity, and the density of the material, respectively. If these coefficients depend on space, one has to use the more general equation  
\[\frac{\partial u}{\partial t}=\;\frac{1}{c\rho}\
abla(k\
abla u), \tag{2}\]  
where it is assumed that the \(c\) and \(\rho\) functions are positive. This equation is now valid for more than one space dimension. According to [PERSON]'s law of cooling, a term \(K(u_{a}-u)\) can describe (free) convective heat transfer [1], where \(u_{a}\) is the ambient temperature (measured in Kelvin), which can be considered as independent from \(u\). On the other hand, according to [PERSON]'s law [2], the heat loss of a unit surface via electromagnetic radiation can be given by a term \(-\sigma u^{4}\). The incoming radiation, the \(Ku_{a}\) heat gain due to the nonzero temperature of the ambient air, as well as other factors such as heat generated by electric currents, can be collected into a so-called heat source term denoted by \(q\). Based on these considerations, the heat conduction Equation (1) can be extended to include heat convection, radiation, and source terms as follows:  
\[\frac{\partial u}{\partial t}=\frac{1}{c\rho}\
abla(k\
abla u)+q-K\cdot u-\sigma \cdot u^{4}, \tag{3}\]  
where the terms \(q\), \(K\cdot u\), and \(\sigma\cdot u^{4}\) should also be non-negative. Several numerical methods, which have been proposed to solve Equations (1) and (2), belong to the wide group of finite difference schemes (FDM) [3; 4; 5]. This can include different kinds of methods of lines [6], where first the space variables are discretized, and then a well-established ordinary differential equation (ODE) solver is employed to solve the obtained ODE system. These methods are usually classified as either explicit or implicit, but occasionally these two approaches are combined [7]. Implicit methods have excellent stability properties; therefore, they are commonly used for these equations [8; 9; 10; 11; 12; 13]. The price for the stability is that a system of algebraic equations must be solved at each time step. This can imply very slow calculations, particularly in cases when the number of space dimensions is more than one, thus one has very large-sized and non-tridiagonal matrices. In these cases, even the most trivial explicit (Euler) time integration can be considerably faster than the implicit one [14]. Moreover, explicit algorithms can be parallelized much more straightforwardly than in the implicit methods. The main obstacle against the rise of the explicit algorithms is that they are typically only conditionally stable, i.e., the solutions can blow up if the time step size is below the so-called Courant-Friedrichs-Lewy (CFL) limit. For example, explicit Runge-Kutta methods can never be unconditionally stable [15] (p. 60). The coefficients \(c\), \(k\), or the diffusion constant can be highly non-uniform in space [16], e.g., when the physical properties have sharp discontinuities at the material boundaries. In these cases, the range of the eigenvalues of the system matrix has several orders of magnitude, the problem is very stiff, the CFL limit can be extremely small, thus the simulation can be unacceptably time-consuming.