is on the order of $10^{-6}$, which would be easily overwhelmed by the numerical truncation errors.

A way to go around the above problems is to enlarge the total energy flux (and thus the superadiabatic gradient) to make the computation feasible. The stellar thermal structure and convective stability are mainly determined by the energy transport mechanism and the superadiabatic gradient. Using the expression of the diffusive radiative flux, we have $$\begin{aligned}
\label{cstr1}
\frac{d T}{d z}&=&-\frac{3\kappa \rho F_{\rm rad}}{4acT^3},\\
\label{cstr2} \frac{d p}{d z}&=&\frac{p}{(\Delta\nabla+\nabla_{\rm
a})T}\frac{d T}{d z}=\rho \vec{g},\end{aligned}$$ where $$\begin{aligned}
\nabla_{\rm a}&=&\nabla_{\rm a}(p,T),\\
\kappa&=&\kappa(p,T),\\
\rho&=&\rho(p,T).\end{aligned}$$ Given the distributions of $F_{\rm rad}$ and $\Delta\nabla$, the EOS and opacity, and the values of $p$, $T$ at one boundary, we can integrate (\[cstr1\], \[cstr2\]) to construct a stratified layer of fluid close to the reference stellar model.

In order to reduce the thermal relaxation time for the intermediate mass star model, the radiative flux $F_{\rm rad}$ (=$F_{\rm b}$ in the radiative region) must be enlarged dramatically (by a factor of $a_{\rm f}\gg 1$). This would make the thermal structure calculated by equations (\[cstr1\], \[cstr2\]) totally different from the reference stellar structure. This problem can be solved by dividing the opacity the same factor $a_{\rm f}$. 