## 4 Stability Analysis of Multivariable Milling Tools
### Third-Order EMHPM for Multivariable Milling Tool  
To approximate the term associated with the delayed terms \(\mathbf{x}_{i}^{\tau_{d}}(T)\) of Equation (29), the interval of the period \(\tau_{T},~{}[t_{0}-\tau_{T},~{}t_{0}]\) is discretized in \(N-1\) intervals that can be equal size. For simplicity, intervals of equal size \(\Delta t=\tau_{T}/(N-1)\) are chosen. Then it is assumed that the function \(\mathbf{x}_{i}^{\tau_{d}}~{}~{}(T)\), which is defined in the interval \([t_{i-d-1},~{}t_{i-d+2}]\), for the third-order EMHPM has the representation of the form:  
\[\begin{array}{l}\mathbf{x}_{i}^{\tau_{d}}(T)=\mathbf{x}_{i-d}(T)\approx \mathbf{x}_{i-d-1}+\big{(}\frac{N-1}{\tau_{T}}\big{)}T\big{(}-\frac{11}{6}( \mathbf{x}_{i-d-1})+3(\mathbf{x}_{i-d})-\frac{3}{2}(\mathbf{x}_{i-  
Then, solving Equation (31) yields  
\[\begin{array}{l}\mathbf{x}_{0}=\mathbf{x}_{i-1}\\ \\ \mathbf{x}_{i1}=T\mathbf{A}_{t}(\mathbf{x}_{i-1})+\sum\limits_{d-\min(d_{n})}^{ \max(d_{n})}\left\{\begin{array}{l}TB_{t}{}^{d}(\mathbf{x}_{i-d-1})+B_{t}{}^ {d}\frac{T^{2}}{2}\Big{(}\frac{N-1}{\tau_{T}}\Big{)}\!\!\left(-\frac{11}{6}( \mathbf{x}_{i-d-1})+3(\mathbf{x}_{i-d})\ -\frac{3}{2}(\mathbf{x}_{i-d+1})+\frac{1}{3}( \mathbf{x}_{i-d+2})\right)\\ +B_{t}{}^{d}\frac{T^{3}}{6}\Big{(}\frac{N-1}{\tau_{T}}\Big{)}^{2}(2(\mathbf{x}_ {i-d-1})-5(\mathbf{x}_{i-d})+4(\mathbf{x}_{i-d+1})-\mathbf{x}_{i-d+2})\\ +B_{t}{}^{d}\frac{T^{3}}{24}\Big{(}\frac{N-1}{\tau_{T}}\Big{)}^{3}(-\mathbf{x} _{i-d-1}+3(\mathbf{x}_{i-d})-3(\mathbf{x}_{i-d+1})+\mathbf{x}_{i-d+2})\\ \end{array}\right\}\\ \\ \mathbf{x}_{2}=\frac{T^{2}}{2}\mathbf{A}_{t}{}^{2}(\mathbf{x}_{i-1})+ \mathbf{A}_{t}\sum\limits_{d-\min(d_{n})}^{\max(d_{n})}\left\{\begin{array}{ l}\frac{T^{2}}{2}B_{t}{}^{d}(\mathbf{x}_{i-d-1})+B_{t}{}^{d}\frac{T^{3}}{6} \Big{(}\frac{N-1}{\tau_{T}}\Big{)}\!\!\left(-\frac{11}{6}(\mathbf{x}_{i-d-1})+ 3(\mathbf{x}_{i-d})\ -\frac{3}{2}(\mathbf{x}_{i-d+1})+\frac{1}{3}(\mathbf{x}_{i-d+2}) \right)\!\!+\\ B_{t}{}^{d}\frac{T^{2}}{24}\Big{(}\frac{N-1}{\tau_{T}}\Big{)}^{2}(2(\mathbf{x}_ {i-d-1})-5(\mathbf{x}_{i-d})+4(\mathbf{x}_{i-d+1})-\mathbf{x}_{i-d+2})\!+\\ B_{t}{}^{d}\frac{T^{5}}{20}\Big{(}\frac{N-1}{\tau_{T}}\Big{)}^{3}(-\mathbf{x}_ {i-d-1}+3(\mathbf{x}_{i-d})-3(\mathbf{x}_{i-d+1})+\mathbf{x}_{i-d+2})\\ \end{array}\right\}\\ \\ \vdots\\ \mathbf{x}_{k}=\frac{1}{k}T^{k}\mathbf{A}_{t}{}^{k}(\mathbf{x}_{i-1})+ \mathbf{A}_{t}{}^{k-1}\sum\limits_{d-\min(d_{n})}^{\max(d_{n})}\left\{ \begin{array}{l}\frac{1}{k}T^{k}\mathbf{B}_{t}{}^{d}(\mathbf{x}_{i-d-1})+ \frac{1}{k+1}T^{k}\mathbf{B}_{t}{}^{d}\Big{(}\frac{N-1}{\tau_{T}}\Big{)}\!\! \left(-\frac{11}{6}(\mathbf{x}_{i-d-1})+3(\mathbf{x}_{i-d})-\frac{3}{2}( \mathbf{x}_{i-d+1})+\right.\\ \frac{1}{3}(\mathbf{x}_{i-d+2})\Big{)}+\frac{1}{(k+2)!}T^{k}{}^{k+2}{}B_{t}{} ^{d}\Big{(}\frac{N-1}{\tau_{T}}\Big{)}^{2}(2(\mathbf{x}_{i-d-1})-5(\mathbf{x}_ {i-d})+4(\mathbf{x}_{i-d+1})-\mathbf{x}_{i-d+2})\!+\\ \left.\frac{1}{(k+3)!}T^{k+3}{}B_{t}{}^{d}\Big{(}\frac{N-1}{\tau_{T}}\Big{)}^{3 }(-\mathbf{x}_{i-d-1}+3(\mathbf{x}_{i-d})-3(\mathbf{x}_{i-d+1})+\mathbf{x}_{i -d+2})\right\end{array}\right\}\\ \\ \mathbf{x}_{k}=\frac{1}{k}T^{k}\mathbf{A}_{t}{}^{k}(\mathbf{x}_{i-1})+ \mathbf{A}_{t}{}^{k-1}\sum\limits_{d-\min(d_{n})}^{\max(d_{n})}\left\{ \begin{array}{l}\frac{1}{k}T^{k}\mathbf{B}_{t}{}^{d}(\mathbf{x}_{i-d-1})+ \frac{1}{k+1}T^{k+1}\mathbf{B}_{t}{}^{d}\Big{(}\frac{N-1}{\tau_{T}}\Big{)}\!\! \left(-\frac{11}{6}(\mathbf{x}_{i-d-1})+3(\mathbf{x}_{i-d-1})+3(\mathbf{x}_{i-d })-\frac{3}{2}(\mathbf{x}_{i-d+1})+\right.\\ \left.\frac{1}{3}(\mathbf{x}_{i-d+2})\right)\!+\!\left(\frac{1}{k+2)!}T^{k+2}{}B_{ t}{}^{d}\Big{(}\frac{N-1}{\tau_{T}}\Big{  
Thus, the stability graphs of Equation (28) were determined by computing the eigenvalues of the transition matrix of Equation (41). The results obtained from the EMHPM were corroborated with the stability lobes in the study of multivariate tools [27].