## 3 Numerical Experiments for the Convection Term  
In Sections 3 and 5, our aim is to mathematically examine the algorithms, thus the units of the quantities are omitted. In this paper, the MATLAB software was used for all numerical calculations. The accuracy is characterized by the usual \(L_{\infty}\) error, which compares the accurate reference value \(u_{i}^{\text{ref}}\) and the result \(u_{i}^{\text{num}}\) obtained by the actual numerical method at the final time \(t^{\text{fin}}\):  
\[Error=\max_{1\leq i\leq N}\Big{|}u_{i}^{\text{ref}}(t^{\text{fin}})-u_{i}^{ \text{num}}(t^{\text{fin}})\Big{|}. \tag{27}\]  
This error is calculated as a function of the time step size \(h\). More concretely, the error is first calculated for a very large \(h\), then this is repeated with time step sizes subsequently decreased by a factor of 2 until small error values are reached.  
Experiment 1--moderately strong convection. We start with a case where the convection coefficient is not very large, \(K=3\), while \(\alpha=1\). The initial and the heat source functions are the following:  
\[u^{0}(x)=5e^{x/4}[2+\sin(2x)]\ \,\ q(x)=10[1+\sin(5x+0.2)]. \tag{28}\]  
The considered domain is \(x_{0}=0\,\ x_{N}=5\,\ t^{0}=1\,\ t^{\text{fin}}=2.2\). The space step size is \(\Delta x=0.0125\). The fixed Dirichlet boundary conditions are simply the initial values at the boundaries \(u^{0}(x_{0})\) and \(u^{0}(x_{N})\). We employed the ode15s solver to gain the reference solution. It is a professionally coded, variable-step, variable-order solver based on the (implicit) numerical differentiation formulas of orders 1 to 5. It is used with a very narrow tolerance (\(10^{-10}\)) to create an accurate reference solution, which is utilized in Equation (27) to compute the maximum error. These errors are presented as a log-log graph in Figure 2. The thin dashed line is exactly proportional to the second power of the time step size, and one can clearly see that the, e.g., ADE and the symmetric mixed version of the LH method have second-order convergence.  
In Figure 3, we also present the initial and final temperatures for the LH method for \(h=9.38\cdot 10^{-3}\). In this case, the maximum error of the \"inside\" implementation is 0.949, while it is 0.926 for the mixed method. One can observe that the numerically obtained functions are very similar to the reference solution, even in the case of these rather large (absolute) errors.  
Figure 3: The temperature \(u\) vs. the \(x\) variable in the case of Experiment 1 for the initial function \(\mathrm{u}^{0}\), the reference solution at \(t^{\mathrm{fin}}\), and the LH method for \(h=2\cdot 10^{-3}\) when the convection term is inside the numerator, and when it is implemented in a mixed way.  
Figure 2: The \(L_{\infty}\) errors as a function of the time step size \(h\) for moderate \(K\) (Experiment 1).  
Experiment 2--very strong convection. Now, a case is considered where the convection coefficient is much larger, \(K=100\), while \(\alpha\) is still one. The initial and the heat source functions are the following:  
\[u^{0}(x)=e^{x}[1+\sin(2x)]\ \,\ \ q(x)=1000[1+\sin(8x+0.2)].\]  
The considered domain and all other circumstances are the same as in the previous experiment, but the space step size is \(\Delta x=0.005\). The error curves are presented in Figure 4.  
In Figure 5, we also present the initial and final temperatures for the LH method for \(h=2\cdot 10^{-3}\), when the maximum error of the \"inside\" implementation is 3.31, while it is 1.35 for the mixed method. The numerical solutions are very similar to the reference solution, even in the case of these rather large errors.  
Figure 4: The \(L_{\infty}\) errors as a function of the time step size \(h\) for large \(K\) (Experiment 2).