## 4 Verification Numerical Examples and Discussion
### Second Example  
In this study, the height of the rectangular cross-section was reduced to \(h=0.1\) m (with \(K_{r}=4.32022\cdot 10^{6}\) Nm obtained from Equation (1)). This value led to the case where \(k_{G}{}^{2}>4\cdot E\text{I}_{z}\cdot k\), consequently meaning that the coefficient \(\lambda_{2}\) became an imaginary value (such cases are possible in the engineering practice, although they are less frequent). This changed not only the numerical values of the coupled GDEs solutions but also their mathematical forms. Therefore, the corresponding solutions were obtained in the form of Equation (6):  
\[v_{GDE,1}(x) =4.0822\cdot 10^{-4}\cdot\text{e}^{-5.5647\cdot x}-7.1761\cdot 10^{-3}\cdot \text{e}^{-1.3372\cdot x}-2.7301\cdot 10^{-7}\cdot\text{e}^{1.3272\cdot x}+2.0557\cdot 10^{-10} \cdot\text{e}^{5.5647\cdot x}-3.3333\cdot 10^{-4}\] \[0\leq x\leq 2\text{ m}\] \[v_{GDE,2}(x) =0.9550\cdot 10^{-4}\cdot\text{e}^{-5.5647\cdot x}-7.2236\cdot 10^{-3} \cdot\text{e}^{-1.3272\cdot x}-3.7995\cdot 10^{-8}\cdot\text{e}^{1.3272\cdot x}+1.9061 \cdot 10^{-20}\cdot\text{e}^{5.5647\cdot x}-3.3333\cdot 10^{-4}\] \[2\leq x\leq 6\text{ m}\] \[v_{GDE,3}(x) =-1.9086\cdot 10^{9}\cdot\text{e}^{-5.5647\cdot x}+0.2974\cdot 10^{-3} \cdot\text{e}^{-1.3372\cdot x}-1.0936\cdot 10^{-9}\cdot\text{e}^{1.3272\cdot x}+1. 1899\cdot 10^{-29}\cdot\text{e}^{5.5647\cdot x}-5.3333\cdot 10^{-4}\] \[6\text{ m}\leq x\leq 10\text{ m}\] \[v_{GDE,4}(x) =2.5682\cdot 10^{19}\cdot\text{e}^{-5.5647\cdot x}-2.1234\cdot\text{e}^{ -1.3272\cdot x}-1.0865\cdot 10^{-9}\cdot\text{e}^{1.3272\cdot x}+5.0959\cdot 10^{-33} \cdot\text{e}^{5.5647\cdot x}-5.3333\cdot 10^{-4}\] \[10\text{ m}\leq x\leq 12\text{ m}\]  
A two cb2 psFE model was also applied here initially, whereas the corresponding stiffness matrices and load vectors were obtained from Equations (25) and (29), respectively. After solving the system of six linear equations, three discrete nodal displacement values as well as rotations were obtained. These values were further coupled with polynomial interpolation functions to additionally obtain transverse displacements between the nodes, as shown in Equations (8) and (9). This problem's results for transverse displacements are presented in Figure 6. The abscissa axis represents the distance from the left end of the beam, and the ordinate axis represents the transverse displacement of the two models. In addition, the red and black lines show the solutions of GDEs and the cb2 psFE model, respectively. It is apparent from Figure 6 that the results of the two FE mesh now exhibit almost perfect agreement at all three nodes and that only some rather small discrepancies are apparent in the vicinities of the cracks. Nevertheless, any discrepancies already completely disappear when four FEs are being implemented, as shown in Figure 6b. Figure 5: Comparison of shear forces for two (**a**) and four (**b**) FE meshes, respectively. Afterwards, several analyses with more finite elements were executed where the displacement curves of all models visually exhibited ideal matching with the GDEs solutions, as no discrepancies were detected. Therefore, just a comparison of displacements at five characteristic points for different finite element meshes is presented in \[0\leq x\leq 2\text{ m}\]0  
Figure 7 (in which the horizontal axis represents the distance from the left end of the beam, and the ordinate axis represents the bending moment) compares the bending moment values obtained by different approaches. As a result, it clearly shows the evolution of bending moment functions obtained from GDEs solutions (red line) as well as interpolated functions obtained by H1 (light blue) and H2 (dark blue) polynomials. It is undoubtedly evident from the figure that the GDEs bending moment functions are dominated by two extremely localized peaks near both beam ends (with intensive changes of curvature within their vicinities).