I have a question regarding the $\textbf{radiation boundary treatment}$ for the 1D heat conduction equation using the $\textbf{implicit}$ finite difference method.
If I use the $\textbf{explicit}$ method, I can do it easily as the equation below. Namely, all the values at the previous time step $i$ is known, so $T^{i+1}$ can easily be calculated, including the internal nodes.
$$\varepsilon\sigma[T_{sur}^4-(T_0^i)^4]+k \frac{T_1^i-T_0^i}{\Delta x}=\rho\frac{\Delta x}{2}C_p\frac{T^{i+1}-T_0^i}{\Delta x}$$
However, if I want to solve it $\textbf{implicitly}$, the equation becomes:
$$\varepsilon\sigma[T_{sur}^4-(T_0^{i+1})^4]+k \frac{T_1^{i+1}-T_0^{i+1}}{\Delta x}=\rho\frac{\Delta x}{2}C_p\frac{T_0^{i+1}-T_0^i}{\Delta x}$$
If I don't have the radiation term $\varepsilon\sigma[T_{sur}^4-(T_0^{i+1})^4]$, then I can easily rewrite the equations, including those for internal nodes, to $[A]\cdot[T^{i+1}]=[B]$. With this, $T^{i+1}$ can be calculated with the matrix inversion.
Nevertheless, if we have the radiative term, then we have a term $(T^{i+1})^4$, which can not be added into the linear equation system $[A]\cdot[T^{i+1}]=[B]$. How could we treat this implicitly then? I do need some help here.