## 2 Materials and Methods
### Mathematical Analysis
#### 2.3.1 Well-Posedness of Model  
Model (4) is a nonlinear system of non-autonomous differential equations. We rely on the numerical methods to obtain a solution. We want to determine whether the solution \(r(t)\) remains bounded for all values of \(t\). Boundedness means the value of \(r(t)\) does not become infinitely large or negative as \(t\) varies. Equation (3) is a linear first-order non-autonomous differential equation, which can be written as  
\[\frac{dr}{dt}=\phi(t)-\xi(t)r(t),\]where \(\phi=\delta(\cos\left(\omega t+\frac{\pi}{6}\right))^{n}\) and \(\xi=\left(\beta-A\left(\frac{R(t)-\overline{R}}{\sigma}\right)\right)\), with initial conditions \(r(0)\geq 0\), \(\forall t\geq 0\). If we solve (5) using an integration factor, then our unique solution becomes  
\[r(t) = r(0)e^{-\left(\int_{0}^{t}\xi(s)ds\right)}+\int_{0}^{t}\phi(u)e^ {-\left(\int_{u}^{t}\xi(s),ds\right)}\,du \tag{6}\] \[\geq 0\text{ for }\ r(0)\geq 0\text{, and }t\geq 0.\]  
Since \(\xi(t)\leq\xi^{M}=2\delta\) and \(\phi(t)\geq\phi^{m}=\beta-A(\frac{R^{M}-\overline{R}}{\sigma})\), the integral in (6) implies  
\[r(t) \leq r(0)e^{-\phi^{m}(t)}+\frac{\xi^{M}}{\phi^{m}}(1-e^{-\phi^{m}(t)} )\leq\max\left\{r(0),\frac{\xi^{M}}{\phi^{m}}\right\},\forall t\geq 0.\]  
Then we obtain  
\[r(t) \leq \max\left\{r(0),\frac{\xi^{M}}{\phi^{m}}\right\}. \tag{7}\]  
Thus, we have the following two lemmas. Lemma 1: _Equation (5) is a linear dynamical system in the following compact interval, which exhibits an attractive nature_  
\[\zeta=\bigg{\{}r\in\mathbb{R}_{+}\mid 0\leq r\leq\frac{\xi^{M}}{\phi^{m}}\bigg{\}}.\]  
Lemma 2: _Equation (3) admits a unique non-negative \(\omega\)-periodic solution \(r^{*}(t)\geq 0\), which is globally attractive in \(\mathbb{R}_{+}\) and is given by_  
\[r^{*}(t)=\frac{e^{-\int_{0}^{t}\phi(s)ds}\int_{0}^{\omega}\xi(v)e^{\int_{v}^{ \omega}\phi(s)ds}dv}{1-e^{-\int_{0}^{\omega}\phi(s)ds}}+e^{-\int_{0}^{t}\phi(s )ds}\bigg{(}\int_{0}^{t}\xi(v)e^{\int_{0}^{\omega}\phi(s)ds}dv\bigg{)}\geq 0. \tag{8}\]  
Once Equation (3) is separately solved, as explained above, we incorporate its solution where appropriate in the remaining equations of the model (4). Since \(r^{*}(t)\geq 0\) is globally attractive in \(\mathbb{R}_{+}\), the system (4) is asymptotically equivalent to the following non-autonomous differential equation system:  
\[ \begin{split}\frac{dP}{dt}&=r^{*}(t)P\bigg{(}1- \frac{P}{K}\bigg{)}-\frac{a}{b+P}LP\\ \frac{dL}{dt}&=L\bigg{[}c\frac{a}{b+P}P-d-\mu\bigg{]},\end{split}  \tag{9}\]  
such that \(P(0)\geq 0\) and \(L(0)\geq 0\). For all state variables \(P(t)\) and \(L(t)\), we should ensure that for all positive times, these state variables must be non-negative for each solution of the system starting with \(P(0)\geq 0\) and \(L(0)\geq 0\). This is equivalent to saying that the positive quadrant \(\mathbb{R}_{2}^{+}=\{(P,L)\in\mathbb{R}_{2}:P\geq 0,L\geq 0\}\) is positively invariant. We denote by \(\mathbb{R}_{2}^{+}\) the non-negative quadrant and by \(int(\mathbb{R}_{2}^{+})\) the positive quadrant. Lemma 3: _(Positivity). All solutions of system (9), starting in \(int(\mathbb{R}_{2}^{+})\), remain positive._  
Proof: From the first equation of system (9), we obtain  
\[\bar{P}(t)=P(0)e^{\int_{0}^{t}\left(r^{*}(t)P\left[1-\frac{P}{K}\right]-\frac {\mu}{P+P}LP\right)ds}\implies P(t)\geq 0. \tag{10}\]From the second equation of the system (4), and using (10), we obtain  
\[\bar{L}(t)=L(0)e^{\left\{\int_{0}^{t}\left(\frac{a}{b+P}-d-\mu\right)ds\right\}} \implies L(t)\geq 0. \tag{11}\]  
If \(P(0)=0\) then \(P(t)=0\) and \(L(t)=0\). That is, if there is no plant, the livestock population will be extinct. If \(L(0)=0\), then \(L(t)=0\) for all \(t>0\), and then the plant population obeys the logistic growth rate law for \(r(t)\). Clearly, exponential functions are always positive. Therefore, if \(P(0)>0\) and \(L(0)>0\), we have the whole dynamics ranging within the first quadrant, that is, \(P(t)>0\) and \(L(t)>0\).