## 2 Problem Statement
Let us note that the introduced dimensionless parameters and functions are expressed through the physical ones, denoted by the containing line above them, with accordance to a rule:  
\[A_{j}=\frac{\bar{A}_{j}}{\sqrt{I_{01}}},\,D_{j}=-\frac{1}{2} \frac{\partial^{2}\bar{k}}{\partial\bar{\omega}^{2}}\Big{|}_{\bar{\omega}_{j}} \frac{Z_{n}}{\tau_{P}^{2}},\omega_{j}=\bar{\omega}_{j}\tau_{P},\,j=1,2,3,\,t= \frac{\bar{t}}{\tau_{P}},\,z=\frac{\bar{z}}{Z_{n}},\] \[\
u_{j1}=\left(\frac{\partial\bar{k}}{\partial\bar{\omega}} \Big{|}_{j\rightarrow\bar{\omega}}-\frac{\partial\bar{k}}{\partial\bar{\omega}} \Big{|}_{\bar{\omega}}\right)\frac{Z_{n}}{\tau_{P}},\,\Delta_{j1}k=\Delta_{j1} \bar{k}Z_{n},\,j=2,\,3 \tag{3}\] \[\gamma=\frac{2\pi\chi^{(2)}\bar{k}\sqrt{I_{01}}}{n^{2}}Z_{n},\]where \(n\) is a refractive index at FF. Parameters \(\bar{k}_{j}\) and \(\Delta_{j1}\bar{k}\) are dimensional wave numbers and phase mismatching, respectively. \(\tau_{P}\) denotes the duration of the incident FW pulse, and \(I_{01}\) denotes its maximum intensity. \(Z_{n}\) is a normalization length equal to 1  mm. For computer simulation, we choose the following values of the physical parameters: \(\tau_{P}\) is changed from 100 fs to 10 ps (meaning the pulses' dispersion actions can be neglected) and the wavelength corresponding to wave with FF is equal to \(\lambda_{1}=1064\)  nm. The angle \(\theta\) of the phase matching for the THG and the second-order susceptibility of selected crystals dependent on the angles \(\theta\) and \(\varphi\) as well as the incident pulse intensity at which the coupling coefficient \(\gamma\) is equal to 1 at the phase matching angle \(\theta\), are shown in Table 1. The angles are measured between the crystal axis in the uni-axial crystal and wave vector of the FW. References to the papers which contain the quadratic susceptibility coefficient are placed after the names of the crystals. Phase matching angles are calculated from the Sellmeier equations, which can be found in [58; 59; 60; 61], respectively. In Table 1, the dimensional phase mismatching between SH and FW is shown. Those values for the potassium beryllium fluoroborate (KBe\({}_{2}\)BO\({}_{3}\)F\({}_{2}\) or KBBF) and rubidium fluoroberyllium borate (RbBe\({}_{2}\)BO\({}_{3}\)F\({}_{2}\) or RBBF) crystals differ insignificantly. We use these values for computer simulation. The computer simulation results obtained for the other three crystals are similar, but they develop at longer distance in comparison with those required using the KBBF and RBBF crystals. It is important to note that at chosen intensities, the cubic nonlinearity only has a weak influence (the corresponding dimensionless coefficient equals a few hundredths), allowing us to neglect it. The set of Equations (1) possess some invariants. The first invariant (the energy preservation law):  
\[I_{1}=\int\limits_{0}^{L_{t}}\Big{(}|A_{1}|^{2}+|A_{2}|^{2}+|A_{3}|^{2}\Big{)}dt \tag{4}\]  
The second invariant:  
\[ \begin{split} I_{2}=\int\limits_{0}^{L_{t}}\Bigg{(}\sum\limits_ {j=1}^{3}p_{j}A_{j}\frac{\partial A_{j}^{*}}{\partial t}\Bigg{)}dt,\,p_{1}=6, \,p_{2}=3,\,p_{3}=2\end{split}  \tag{5}\]  
The third invariant (Hamiltonian):  
\[ \begin{split} I_{3}=\int\limits_{0}^{L_{t}}\Big{(}-\sum\limits_ {j=1}^{3}p_{j}D_{j}\Big{|}\frac{\partial A_{j}}{\partial t}\Big{|}^{2}-p_{2} \
u_{21}\,\text{Im}\Big{(}A_{2}^{*}\frac{\partial A_{2}}{\partial t}\Big{)}- p_{3}\
u_{31}\,\text{Im}\Big{(}A_{3}^{*}\frac{\partial A_{3}}{\partial t}\Big{)}\\ +6\gamma\,\text{Re}\big{(}2A_{1}A_{2}A_{3}^{*}e^{i(\Delta_{31}k- \Delta_{2}k)z}+A_{1}^{2}A_{2}^{*}e^{i\Delta_{2}k+z}\big{)}+3\Delta_{21}k|A_{2} |^{2}+2\Delta_{31}k|A_{3}|^{2}\Big{)}dt,\end{split}  \tag{6}\]