## 3 The Procedure of the Numerical Solution  
### The Spatial and Temporal Discretization  
We consider the case of heat conduction because the concepts of the following discretization process are better established in this case. Let us discretize the time variable uniformly, which means \(t\in\left[t^{0},\;\;t^{\mathrm{fin}}\right]\), and  
\[t^{n}=t^{0}+nh\;,\;\;n=1,...,\;\;T\;,\;\;hT=t^{\mathrm{fin}}-t^{0}\]  
On the interval \(x\in\left[x_{0}\;,\;x_{N}=x_{0}+L\right]\subset\mathbb{R}\), we construct an equidistant spatial grid:  
\[x_{j}=x_{0}+j\Delta x\;,\;\;j=0,...,\;\;N\;,\;\;N\Delta x=L\]  
The parameters \(t^{0}\), \(x_{0}\), etc., will be given when the concrete examples are presented.  
If the physical properties of the heat-conducting media depend on space, the following PDE is used:  
\[c(x)\rho(x)\frac{\partial u}{\partial t}=\ \ \frac{\partial}{\partial x}\left(k(x) \frac{\partial u}{\partial x}\right) \tag{7}\]  
For simplicity, we consider \(c(x,t)\equiv 1\) and \(\rho(x,t)\equiv 1\) in this work, and all the space and time dependence of the diffusivity will be incorporated into the conductivity \(k(x,t)\). When the outer differentiation with respect to \(x\) is executed at the right-hand side of (7), the term \(k(x)\) must also be differentiated, which implies that an extra drift term with \(\partial u/\partial x\) appears, such as in Equation (2). We avoid this by discretizing the function \(k\) and simultaneously \(\partial u/\partial x\) in Equation (7). The usual central difference formula is employed to obtain  
\[\frac{\partial u}{\partial t}\bigg{|}_{x_{i},t^{n}}=\frac{1}{\Delta x}\left[k \left(x_{i}+\frac{\Delta x}{2},t^{n}\right)\frac{u(x_{i}+\Delta x)-u(x_{i})}{ \Delta x}+k\left(x_{i}-\frac{\Delta x}{2},t^{n}\right)\frac{u(x_{i}-\Delta x)- u(x_{i})}{\Delta x}\right]\]  
We now move from node variables to cell variables. It means that \(u_{i}\) is the approximation of the average temperature of cell \(i\), by its value at the cell center. Furthermore, \(k_{i,j+1}\) is the heat conductivity between cell \(i\) and its (right) neighbor, estimated by its value at the border of the cells. Now, the previous formula will have the form  
\[\frac{du_{i}}{dt}=\frac{1}{\Delta x}\left(k_{i,i+1}^{n}\frac{u_{i+1}-u_{i}}{ \Delta x}+k_{i,i-1}^{n}\frac{u_{i-1}-u_{i}}{\Delta x}\right)\]  
Because we are in one spatial dimension, we consider the cross-section area of the system as unity. Now, the heat capacity of the cell is the same as the volume and can be given as \(C_{i}=V_{i}=\Delta x\). The thermal resistances at the \(n\)-th time level are calculated as follows:  
\[R_{i,\ i+1}^{n}=\frac{\Delta x}{k_{i,i+1}^{n}}=\ \frac{\Delta x}{D\left(x_{i,i+1} /\sqrt{t^{n}}\right)}m,\ i=1,\ldots,N-1\]  
Now, we have the equation for the time derivative of each cell variable:  
\[\frac{du_{i}}{dt}=\frac{u_{i-1}-u_{i}}{R_{i,\ i-1}^{n}C_{i}}+\frac{u_{i+1}-u_ {i}}{R_{i,\ i+1}^{n}C_{i}} \tag{8}\]  
which can be written into a matrix form  
\[\frac{d\overrightarrow{u}}{dt}=M\overrightarrow{u} \tag{9}\]  
where the system matrix \(M\) is \(N\times N\) dimensional and depends on the time variable. More details about this way of discretization (for the case of the time-independent case but more space dimensions) can be found, e.g., in [33].  
The numerical methods examined in this work can be used in the case of different types of boundary conditions, e.g., Dirichlet, Neumann, and periodic. Because we reproduce the above given analytical solutions, we always use Dirichlet boundary conditions, but the concrete details will be given later.