## 2 The Examined Numerical Methods
### Implementations of the Radiation Term  
In Sections 3 and 4, it will turn out that the mixed treatment with \(p=\
icefrac{{1}}{{2}}\) is the most effective; thus, when defining different treatments of the radiation term, the convection term will be taken into account that way. The radiation term will be implemented in similar ways as the convection term, but now one has more possibilities. 1. Explicit treatment:  
\[u_{i}^{\sfrac{1}{2}}=\frac{u_{i}^{0}+\
icefrac{{r}}{{2}}\big{(}u_{i-1}^{0}+u _{i+1}^{0}\big{)}+\
icefrac{{h}}{{2}}\cdot q_{i}-Khu_{i}^{0}/2}{1+r+Kh/4}- \sigma h\big{(}u_{i}^{0}\big{)}^{4}/2. \tag{16}\]  
2. Quasi-exact treatment: The analytical solution of the ODE  
\[\frac{du}{dt}=-\sigma\cdot u^{4},\]  
is  
\[u(t)=\Big{(}(u_{0})^{-3}+3\sigma t\Big{)}^{-\frac{1}{3}}. \tag{17}\]  
This means that we have the following two sub-stages:  
\[ \begin{split}\text{Sub}-\text{stage}\;1:&\;u_{i}^{ \sfrac{1}{2},\text{ temp}}=\frac{u_{i}^{0}+\
icefrac{{r}}{{2}}\big{(}u_{i-1}^{0}+u_{i+1}^{0} \big{)}+\
icefrac{{h}}{{2}}\cdot q_{i}-Khu_{i}^{0}/2}{1+r+Kh/4},\\ \text{Sub}-\text{stage}\;2:&\;u_{i}^{\sfrac{1}{2}}= \bigg{(}\big{(}u_{i}^{\sfrac{1}{2},\text{ temp}}\big{)}^{-3}+3\sigma h\bigg{)}^{-\frac{1}{3}}.\end{split}  \tag{18}\]  
3. Pseudo-implicit treatment: Equation (4) is now modified as  
\[u_{i}^{n+1}=u_{i}^{n}+r\Big{[}u_{i-1}^{n}-2u_{i}^{n}+u_{i+1}^{n}-2\theta u_{i} ^{n}-2(1-\theta)u_{i}^{n+1}\Big{]}-Khu_{i}^{n+1}/2-Khu_{i}^{n}/2-\sigma hu_{i} ^{n+1}(u_{i}^{n})^{3} \tag{19}\]  
This yields  
\[u_{i}^{\sfrac{1}{2}}=\frac{u_{i}^{0}+\
icefrac{{r}}{{2}}\big{(}u_{i-1}^{0}+u_ {i+1}^{0}\big{)}+\
icefrac{{h}}{{2}}\cdot q_{i}-Khu_{i}^{0}/2}{1+r+Kh/4+\sigma h \big{(}u_{i}^{0}\big{)}^{3}/2} \tag{20}\]4. 'Inside' treatment: The last term of Equation (19) is now written as \(-\sigma h\big{(}u_{i}^{n}\big{)}^{4}\), which yields  
\[u_{i}^{\
icefrac{{1}}{{2}}}=\frac{u_{i}^{0}+\
icefrac{{r}}{{2}}\big{(}u_{i-1}^{0}+u _{i+1}^{0}\big{)}+\
icefrac{{h}}{{2}}\cdot q_{i}-Khu_{i}^{0}/2-\sigma h\big{(}u_{ i}^{0}\big{)}^{4}/2}{1+r+Kh/4}. \tag{21}\]  
5. Mixed treatment with equal share of the pseudo-implicit and inside treatments. The last term of Equation (19) is the average of the previous two cases, i.e., it is \(-\sigma h\big{(}u_{i}^{n}\big{)}^{4}/2-\sigma h\mu_{i}^{n+1}\big{(}u_{i}^{n} \big{)}^{3}/2\), which yields  
\[u_{i}^{\
icefrac{{1}}{{2}}}=\frac{u_{i}^{0}+\
icefrac{{r}}{{2}}\big{(}u_{i-1 }^{0}+u_{i+1}^{0}\big{)}+\
icefrac{{h}}{{2}}\cdot q_{i}-Khu_{i}^{0}/2-\sigma h \big{(}u_{i}^{0}\big{)}^{4}/4}{1+r+Kh/4+\sigma h\big{(}u_{i}^{0}\big{)}^{3}/4} \tag{22}\]  
Now, we turn our attention to the nonstandard or nonlocal treatments of the radiation term. Aiming to avoid symmetry breakings and an extensive increase in the running times, we try three different possibilities. Since the pseudo-implicit version (20) is the most successful among the treatments presented so far, we modify this version, mostly by changing one or two of the \(u_{i}\)-s in the \(\left(u_{i}\right)^{3}\) product in the denominator of, e.g., Equation (20) in the following three ways. 6. Product treatment (denoted by LH PI NL prod): Instead of \(\left(u_{i}^{0}\right)^{3}\) and \(\left(u_{i}^{n}\right)^{3}\), we write \(u_{i-1}^{0}u_{i}^{0}u_{i+1}^{0}\) and \(u_{i-1}^{n}u_{i}^{n}u_{i+1}^{n}\), respectively. 7. Average treatment (denoted by LH PI NL av): Instead of \(\left(u_{i}^{0}\right)^{3}\) and \(\left(u_{i}^{n}\right)^{3}\), we write \(\frac{u_{i-1}^{0}+u_{i+1}^{0}}{2}\big{(}u_{i}^{0}\big{)}^{2}\) and \(\frac{u_{i-1}^{n}+u_{i+1}^{n}}{2}\big{(}u_{i}^{n}\big{)}^{2}\), respectively. 8. In the case of the time-average treatment (denoted by LH PI NL time), there are two sub-stages. First, we calculate the effect of the diffusion and the source terms by (5)-(7) to obtain a temporary value \(u_{i}^{\text{temp}}\).