## Appendix C Numerical Scheme  
The numerical integration scheme employed the backward differentiation formula (BDF) of the second order with an automatic step-size selection. In order to understand the numerical scheme, consider an ordinary differential equation  
\[\frac{\partial x}{\partial t}=f(t,x),\text{ for }0<t<T, \tag{10}\]  
given initial point \(x(0)=x_{0}\). Let us use [PERSON]'s notation and denote \(\vec{z}(t)\) the vector of time-derivatives with a step size \(h>0\), that is,  
\[\vec{z}=\left(x,\,h\dot{x},\,\frac{h^{2}}{2}\vec{x}\right). \tag{11}\]  
The first two components of \(\vec{z}_{0}\) are \(x_{0}\) and \(f(0,x_{0})\), while the second derivative should be set to zero as we assumed the transient started from a steady state.  
Using [PERSON]'s notation, the multistep BDF can be formulated as a predictor-corrector scheme. The predictor \(\vec{z}_{h}^{p}\) of [PERSON]'s vector at the time point \(t+h\) is given by  
\[\vec{z}_{h}^{p}=\begin{bmatrix}1&1&1\\ 0&1&2\\ 0&0&1\end{bmatrix}\vec{z}_{0}. \tag{12}\]This equation can be extended to an arbitrary order considering the matrix on the right-hand side is the Pascal matrix. The correction step applied to the vector \(\vec{z}^{P}\) is  
\[\vec{z}_{h}=\vec{z}_{h}^{P}+\vec{l}\left(x_{h}-x_{h}^{P}\right). \tag{100}\]  
Here, coefficients \(\vec{l}\) depend on the numerical scheme order. For the second order, \(l_{0}=0\), \(l_{1}=3/2\) and \(l_{2}=1/2\). The solution \(x(t+h)\) is calculated as follows  
\[x_{h}=x_{h}^{P}+hl_{1}^{-1}\left(f(t+h,x_{h})-\dot{x}_{h}^{P}\right). \tag{101}\]  
The local truncation error \(e_{h}\) of solution \(x_{h}\) at the time point \(t+h\) relates to the difference between the predictor and corrector vectors \(\vec{z}_{h}^{P}\) and \(\vec{z}_{h}\). Specifically, for the \(n\)th order scheme, the local error is given by  
\[e_{h}=l_{n}\left(x_{h}-x_{h}^{P}\right)+O(h^{n+1}). \tag{102}\]  
The time step is selected in such a way that the relative local error is less than one. Given absolute and relative tolerances \(A\) and \(R\), the step size \(h^{\prime}\) is  
\[h^{\prime}=C\left(\frac{A+Rx_{h}}{x_{h}-x_{h}^{P}}\right)^{\frac{1}{n+1}}h, \tag{103}\]  
where the safety coefficient \(C=1.2\).  
The extension to the system of ODEs is straightforward. [PERSON]'s vector \(\vec{z}\) is replaced by a matrix whose columns correspond to one equation of the ODE system, the division in Equation (103) is assumed component-wise and the brackets in that equation are replaced by a norm. The other equations of the integration scheme do not change.  
The fixed-point algorithm iterates the solution of nonlinear system Equation (101) to achieve the ultimate convergence for the given floating-point number precision. Specifically, the iteration stops when all the components of solution vector \(\vec{x}_{h}\) become periodical. It is clear that for such a stopping criterion, the higher the floating-point precision used, the more iterations have to be done until convergence. However, the convergence rate of the fixed-point iteration algorithm for the point-kinetic model is rather high; therefore, even for the double precision, the number of iterations is less than 20.