## 2 Semi-Analytical Solution of the Duffing Model  
The forced and damped Duffing oscillator expression is as follows:  
\[\ddot{x}+\delta\dot{x}-\alpha x+\beta x^{3}=Q_{0}\cos\Omega t \tag{1}\]where constant \(\alpha>0\), \(\beta>0\) is the stiffness coefficient of soft spring. Here, we analyze the bifurcations of the Holmes-Duffing system, which is a typical Duffing oscillator. Let \(\alpha\), \(\beta=1\). \[\ddot{x}+\delta\dot{x}-x+x^{3}=Q_{0}\cos\Omega t \tag{2}\]  
The state equation of Equation (2) in state space is  
\[ \begin{split}\dot{x}&=y\\ \dot{y}&=Q_{0}\cos\Omega t-(\delta\dot{x}-\alpha x+ \beta x^{3}).\end{split}  \tag{3}\]  
From the work of [PERSON] and [PERSON] [24; 25], the discretization of Equation (3) can be achieved by applying a midpoint scheme during the time interval \([t_{k-1},t_{k}]\), resulting in the following implicit expression:  
\[ \begin{split} x_{k}&=x_{k-1}+\frac{1}{2}h\big{(}y_{k -1}+y_{k}\big{)},\\ y_{k}&=y_{k-1}+h[Q_{0}\cos\Omega(t_{k-1}+\frac{1}{ 2}h)-\frac{1}{2}\delta(y_{k-1}+y_{k})\\ &+\frac{1}{2}\alpha(x_{k-1}+x_{k})-\frac{1}{8}\beta(x_{k-1}+x_{k })^{3}].\end{split}  \tag{4}\]  
where \(t_{0}\) represents the initial time, and \(h=t_{k}-t_{k-1}\) is the discrete time step. The structure of the implicit mapping can be depicted as follows:  
\[ \begin{split}& P_{1}:(x_{0},y_{0})\rightarrow(x_{1},y_{1}) \Longrightarrow(x_{1},y_{1})=P_{1}(x_{0},y_{0});\\ & P_{2}:(x_{1},y_{1})\rightarrow(x_{2},y_{2})\Longrightarrow(x _{2},y_{2})=P_{2}(x_{1},y_{1});\\ &\vdots\\ & P_{N-1}:(x_{N-2},y_{N-2})\rightarrow(x_{N-1},y_{N-1})\\ &\implies(x_{N-1},y_{N-1})=P_{N}(x_{N-2},y_{N-2});\\ & P_{N}:(x_{N-1},y_{N-1})\rightarrow(x_{N},y_{N})\\ &\implies(x_{N},y_{N})=P_{N}(x_{N},y_{N});\\ &\vdots\\ & P_{mN}:(x_{mN-1},y_{mN-1})\rightarrow(x_{mN},y_{mN})\\ &\implies(x_{mN},y_{mN})=P_{mN}(x_{mN-1},y_{mN-1});\end{split}  \tag{5}\]  
The motion of period-1 can be discretized by \(N\) points, while period-m motions require \(mN\) points. By applying corresponding periodic conditions, we obtain the following expression:  
\[(x_{mN},y_{mN})=(x_{0},y_{0}) \tag{6}\]  
The \(2(mN+1)\) equations in Equations (5) and (6) are used to explore periodic motions. By calculating the numerical outcomes of such periodic movements, one can subsequently examine the stability and bifurcations by analyzing the eigenvalues of the related Jacobian matrix. Let \(\mathbf{x}_{k}=(x_{k},y_{k})^{T}\) represent a value from the semi-analytical solution, and \(\mathbf{f}_{k}=(f_{1,k},f_{2,k})^{T}\) be the homogenized version of discretization, where \(k=0,1,2,...,N\). In the small perturbations around period-m motion \(\mathbf{x}_{k}^{*}\), the linear equation is derived as  
\[\left.\frac{\partial\mathbf{f}_{k}}{\partial\mathbf{x}_{k-1}}\right|_{( \mathbf{x}_{k-1}^{*}\mathbf{x}_{k}^{*})}\Delta\mathbf{x}_{k-1}+\left.\frac{ \partial\mathbf{f}_{k}}{\partial\mathbf{x}_{k}}\right|_{(\mathbf{x}_{k-1}^{*} \mathbf{x}_{k}^{*})}\Delta\mathbf{x}_{k}=\mathbf{0} \tag{7}\]So, we can have  
\[ \begin{split} DP_{k}&=\frac{\partial\mathbf{x}_{k}}{ \partial\mathbf{x}_{k-1}}=-\left(\frac{\partial\mathbf{f}_{k}}{\partial \mathbf{x}_{k}}\right)^{-1}\frac{\partial\mathbf{f}_{k}}{\partial\mathbf{x}_{ k-1}}\Bigg{|}_{(\mathbf{x}_{k-1}^{\star}\mathbf{x}_{k}^{\star})}\\ &=\left[\begin{array}{cc}\frac{\partial\mathbf{x}_{k}}{ \partial\mathbf{x}_{k-1}}&\frac{\partial\mathbf{x}_{k}}{\partial\mathbf{y}_{ k-1}}\\ \frac{\partial\mathbf{y}_{k}}{\partial\mathbf{x}_{k-1}}&\frac{\partial \mathbf{y}_{k}}{\partial\mathbf{y}_{k-1}}\end{array}\right]\Bigg{|}_{( \mathbf{x}_{k-1}^{\star},\mathbf{x}_{k}^{\star})}\end{split}  \tag{8}\]  
Meanwhile,  
\[\Delta\mathbf{x}_{N}=DP\Delta x_{0}=DP_{N}\cdot DP_{N-1}\cdots DP_{2}\cdot DP_{ 1}\Delta\mathbf{x}_{0} \tag{9}\]  
To analyze the periodic motion, the eigenvalues of the Jacobian matrix, which represent the discretized system, are computed.