## 4 Numerical Method for the ODEs
### Runge-Kutta Fourth-Order Method (Classical)  
Runge-Kutta methods propagate a solution over an interval by combining the information from several Euler-style steps and then using the information obtained to match a Taylor series expansion up to some higher order. The common opinion of Runge-Kutta is that it is the best function to apply as a first try for most problems or for trivial problems where computational efficiency is of no concern and succeeds virtually very often. However, it is not usually the fastest method, except when moderate accuracy 1e-5 is required [30].  
The classical Runge-Kutta method is based on the Simpson quadrature rule and uses four explicit stages, hence four function evaluations per step, to achieve \(O(dt^{4})\) accuracy. Local truncation error is fourth order if \(F(t)\) has five continuous derivatives. In general, the classical fourth-order Runge-Kutta method requires four evaluations \(K\) of the right-hand side of the differential equation per time step [25; 26; 27].  
The Runge-Kutta classical fixed time step size method for the discussed set of ODE balance equations was written as  
\[K_{1}\Big{(}T_{f}\Big{)}=dt\;dT_{f}\Big{(}T_{f},T_{m}\Big{)} \tag{28}\] \[K_{1}(T_{m})=dt\;dT_{m}\Big{(}T_{f},T_{m},\;T_{a}\Big{)}\;,\] (29) \[K_{1}(T_{a})=dt\;dT_{a}(T_{m},T_{a},\;a)\] (30) \[K_{1}(a)=dt\;da(T_{a},a) \tag{31}\]  
\[K_{2}\Big{(}T_{f}\Big{)}=dt\;dT_{f}\Big{(}T_{f}+0.5K_{1}\Big{(}T_{f}\Big{)},T_ {m}+0.5K_{1}(T_{m})\Big{)} \tag{32}\] \[K_{2}(T_{m})=dt\;dT_{m}\Big{(}T_{f}+0.5K_{1}\Big{(}T_{f}\Big{)},T _{m}+0.5K_{1}(T_{m}),\;T_{a}+0.5K_{1}(T_{a})\Big{)}\;,\] (33) \[K_{2}(T_{a})=dt\;dT_{a}(T_{m}+0.5K_{1}(T_{m}),T_{a}+0.5K_{1}(T_{ a}),\;a+0.5K_{1}(a))\] (34) \[K_{2}(a)=dt\;da(T_{a}+0.5K_{1}(T_{a}),a+0.5K_{1}(a)) \tag{35}\]  
\[K_{3}\Big{(}T_{f}\Big{)}=dt\;dT_{f}\Big{(}T_{f}+0.5K_{2}\Big{(}T_{f}\Big{)},T _{m}+0.5K_{2}(T_{m})\Big{)} \tag{36}\] \[K_{3}(T_{m})=dt\;dT_{m}\Big{(}T_{f}+0.5K_{2}\Big{(}T_{f}\Big{)},T _{m}+0.5K_{2}(T_{m}),\;T_{a}+0.5K_{2}(T_{a})\Big{)}\;,\] (37) \[K_{3}(T_{a})=dt\;dT_{a}(T_{m}+0.5K_{2}(T_{m}),T_{a}+0.5K_{2}(T_{ a}),\;a+0.5K_{2}(a))\] (38) \[K_{3}(a)=dt\;da(T_{a}+0.5K_{2}(T_{a}),a+0.5K_{2}(a)) \tag{39}\]  
\[K_{4}\Big{(}T_{f}\Big{)}=dt\;dT_{f}\Big{(}T_{f}+K_{3}\Big{(}T_{f}\Big{)},T_ {m}+K_{3}(T_{m})\Big{)} \tag{40}\] \[K_{4}(T_{m})=dt\;dT_{m}\Big{(}T_{f}+K_{3}\Big{(}T_{f}\Big{)},T_ {m}+K_{3}(T_{m}),\;T_{a}+K_{3}(T_{a})\Big{)}\;,\] (41) \[K_{4}(T_{a})=dt\;dT_{a}(T_{m}+K_{3}(T_{m}),T_{a}+K_{3}(T_{a}),\; a+K_{3}(a))\] (42) \[K_{4}(a)=dt\;da(T_{a}+K_{3}(T_{a}),a+K_{3}(a)) \tag{43}\]  
\[\left\{\begin{array}{l}T_{f(n+1)}=T_{f}+\Big{(}K_{2}\Big{(}T_{f}\Big{)}+2K_ {2}\Big{(}T_{f}\Big{)}+2K_{3}\Big{(}T_{f}\Big{)}+K_{4}\Big{(}T_{f}\Big{)}\Big{)} /6\\ T_{m(n+1)}=T_{m}+(K_{2}(T_{m})+2K_{2}(T_{m})+2K_{3}(T_{m})+K_{4}(T_{m}))/6\,, \end{array}\right. \tag{44}\] \[T_{a(n+1)}=T_{a}+(K_{2}(T_{a})+2K_{2}(T_{a})+2K_{3}(T_{a})+K_{4} (T_{a}))/6\] (45) \[a(n+1)=a+(K_{2}(a)+2K_{2}(a)+2K_{3}(a)+K_{4}(a))/6 \tag{46}\]  
where \(dt\) is the increase of the independent variable \(t\), and \(dT_{f},\;dT_{m},\;dT_{a},\;da\) are time derivatives of searched functions implicitly dependent on time and where \(n\) is the time point of the known values of the function and \(n+1\) is the next time point. The resulting equations were implemented manually as a code and solved in MATLAB(r) platform.