## 2 Analytical Model and the Development of the Mathematica Package
The model treats the cross-section as in-plane rigid and excludes shear effects on its median surface [32; 44]. \[ \begin{array}{c}\Pi=\frac{1}{2}\int\limits_{0}^{L}EI_{y}(z)\Big{[}\frac{d^{2 }u(z)}{dz^{2}}\Big{]}^{2}dz+\frac{1}{2}\int\limits_{0}^{L}GJ(z)\Big{[}\frac{d \phi(z)}{dz}\Big{]}^{2}dz\\ +\frac{1}{2}\int\limits_{0}^{L}E\bigg{[}C_{w}(z)\Big{[}\frac{d^{2}\phi(z)}{dz^ {2}}\Big{]}^{2}+2\Big{[}\frac{d^{2}C_{w}(z)}{dz^{2}}\Big{]}\Big{[}\frac{d\phi (z)}{dz}\Big{]}^{2}+2\Big{[}\frac{d\mathcal{C}_{w}(z)}{dz}\Big{]}\Big{[}\frac {d\phi(z)}{dz}\Big{]}\Big{]}dz\\ +\frac{1}{2}\int\limits_{0}^{L}M_{x}(z)\bigg{[}2\phi(z)\Big{[}\frac{d^{2}u(z) }{dz^{2}}\Big{]}+\beta_{x}(z)\Big{[}\frac{d\phi(z)}{dz}\Big{]}^{2}\bigg{]}dz+ \frac{1}{2}\sum\limits PH_{p}\phi_{p}^{2}+\frac{1}{2}\int\limits_{0}^{L}qH_{q} \phi(z)^{2}dz\end{array}  \tag{16}\]  
The Ritz method relies on the principle that the extra strain energy accumulated during lateral-torsional buckling equals the additional work done by the applied loads. This energy-based approach necessitates an assumed buckled shape that satisfies the kinematic boundary conditions and closely resembles the beam's actual deformed configuration. Consequently, the LTB load is then obtained by substituting the predefined buckling shape into the total potential energy expression. Kinematic boundary conditions, which impose geometric restrictions on displacements and rotations at support points, play a crucial role in defining the admissible shape function. For a simply supported beam, the restraint conditions at the supports (\(z=0\) and \(z=L\)) are \(u=\Phi=0\) and \(d^{2}u/dz^{2}=d^{2}\Phi/dz^{2}=0\). The buckled shape of a web-tapered mono-symmetric I-section beam is described by the lateral displacement of the shear center (\(u\)) and the twist angle of the cross-section (\(\Phi\)). Taking into account the corresponding boundary conditions, the buckled shape functions are constructed as shown in Equation (17), where \(A_{i}\) and \(B_{i}\) represent the amplitudes of the mode shape functions. \[ \begin{array}{c}u(z)=\sum\limits_{i=1}^{n}A_{i}\ sin\Big{[}\frac{i_{i2}}{L} \Big{]}\\ \phi(z)=\sum\limits_{i=1}^{n}B_{i}\ sin\Big{[}\frac{i_{i2}}{L}\Big{]}\end{array}  \tag{17}\]  
By substituting the buckled shape functions defined in Equation (17) into the total potential energy expression in Equation (16), one obtains:  
\[\Pi=U+V=f(A_{1},\ A_{2},\ A_{3}\ldots,\ A_{n},\ B_{1},\ B_{2},\ B_{3}\ldots,\ B_{n}) \tag{18}\]In Equation (18),\(f\) is a function with \(2n\) variables, namely the coefficients \(A_{i}\) and \(B_{i}\). At the onset of lateral-torsional buckling, the total potential energy reaches a stationary value, leading to the following necessary conditions:  
\[ \begin{array}{c}\frac{\partial\Pi}{\partial A_{i}}=0\ \ for\ i=1,\ 2\ldots\ n\\ \frac{\partial\Pi}{\partial B_{i}}=0\ \ for\ i=1,\ 2\ldots\ n\end{array}  \tag{19}\]  
These conditions yield a set of \(2n\) homogeneous linear equations, which can be conveniently expressed in matrix notation:  
\[\left[K\right]_{2n\times 2n}\left\{d\right\}_{2n\times 1}=\left\{0\right\}_{2n \times 1} \tag{20}\]  
where \(K\) is the stiffness-like coefficient matrix, and \(\left\{d\right\}=\left\{A_{1}\text{, }A_{2}\text{, }A_{3}\text{, }\ldots\text{,}A_{n}\text{, }B_{1}\text{, }B_{2}\text{, }B_{3}\text{, }\ldots\text{,}B_{n}\right]^{T}\) contains the unknown amplitudes. A non-trivial solution to Equation (20) exists only when the determinant of \(K\) vanishes:  
\[\text{det}\left[K\right]_{2n\times 2n}=0 \tag{21}\]  
Solving Equation (21) leads to a polynomial equation of degree 2n, whose smallest positive root corresponds to the critical elastic lateral-torsional buckling load or moment. Although this analytical procedure could not establish a closed-form simplified equation, it can be easily programmed by any mathematical software.