**Proof.** See [7].

**Theorem 7.** *The disease-free equilibrium point Pf* = 0, 0, <sup>Λ</sup> *μ*0 *of the scheme (34) is locally asymptotically stable if* R<sup>0</sup> < 1*.*

**Proof.** Calculating the eigenvalues of the Jacobian matrix of the system (34) at the diseasefree point, we obtain the following characteristic polynomial

$$
\left(\frac{1}{1+\varrho(h)\mu\_0} - \lambda\right) \left| \begin{array}{c} \frac{1}{1+\varrho(h)\delta\_I} - \lambda \\ \frac{\mu\_0 \hbar \mathcal{N} \delta\_I}{\mu\_0 \hbar \mathcal{N} \delta\_I} \end{array} \frac{\varrho(h)\beta \varepsilon^{-\delta\_I \Delta} \Lambda}{\mu\_0 (1+\varrho(h)\delta\_I)} \right|\_{\mu\_0} = 0.
$$

Thus, the first eigenvalue is

$$
\lambda\_1 = \frac{1}{1 + \varphi(h)\mu\_0} < 1.
$$

The other two eigenvalues, are the roots of quadratic polynomial

$$
\lambda^2 - \left(\frac{1}{p} + \frac{\mu\_0}{q}\right)\lambda + \frac{\mu\_0 - \varrho(h)^2 \delta\_I N \beta \Lambda e^{-\delta\_{I\_E}\Lambda}}{\mu\_0 pq} = 0,\tag{36}
$$

where *<sup>p</sup>* <sup>=</sup> <sup>1</sup> <sup>+</sup> *<sup>ϕ</sup>*(*h*)*δ<sup>I</sup>* <sup>&</sup>gt; 0 and *<sup>q</sup>* <sup>=</sup> *<sup>μ</sup>*<sup>0</sup> <sup>+</sup> *<sup>μ</sup>*0*ϕ*(*h*)*<sup>C</sup>* <sup>+</sup> *<sup>ϕ</sup>*(*h*)*β*Λ. Next, let *<sup>a</sup>*<sup>1</sup> <sup>=</sup> <sup>1</sup> *p* <sup>+</sup> *<sup>μ</sup>*<sup>0</sup> *q* and *<sup>a</sup>*<sup>2</sup> <sup>=</sup> *<sup>μ</sup>*<sup>0</sup> <sup>−</sup> *<sup>ϕ</sup>*(*h*)2*δINβ*Λ*<sup>e</sup>* −*δIE*<sup>Δ</sup> *<sup>μ</sup>*<sup>0</sup> *pq* . We have the following affirmations. 1. If 1 <sup>&</sup>gt; <sup>R</sup>0, it follows that *<sup>ϕ</sup>*(*h*)2*δI*(*Cμ*<sup>0</sup> <sup>+</sup> *<sup>β</sup>*Λ) <sup>&</sup>gt; *<sup>ϕ</sup>*(*h*)2*δINβ*Λ*<sup>e</sup>* <sup>−</sup>*δIE*Δ. Thus, *δI*(*μ*<sup>0</sup> + *μ*0*hC* + *hβ*Λ) > *μ*0*δ<sup>I</sup>* + *ϕ*(*h*)*δINβ*Λ*e* <sup>−</sup>*δIE*<sup>Δ</sup> ⇐⇒ (1 + *ϕ*(*h*)*δI*)*q* > *q* + *ϕ*(*h*)*μ*0*δ<sup>I</sup>* + *ϕ*(*h*)2*δINβ*Λ*e* <sup>−</sup>*δIE*<sup>Δ</sup> ⇐⇒ *pq* + *μ*<sup>0</sup> > *q* + *ϕ*(*h*)*μ*0*δ<sup>I</sup>* + *ϕ*(*h*)2*δINβ*Λ*e* <sup>−</sup>*δIE*<sup>Δ</sup> <sup>+</sup> *<sup>μ</sup>*<sup>0</sup> ⇐⇒ *pq* + *μ*<sup>0</sup> > *q* + *ϕ*(*h*)2*δINβ*Λ*e* <sup>−</sup>*δIE*<sup>Δ</sup> <sup>+</sup> *<sup>p</sup>μ*<sup>0</sup> ⇐⇒ 1 −  1 *p* <sup>+</sup> *<sup>μ</sup>*<sup>0</sup> *q* <sup>+</sup> *<sup>μ</sup>*<sup>0</sup> <sup>−</sup> *<sup>h</sup>*2*δINβ*Λ*<sup>e</sup>* −*δIE*<sup>Δ</sup> *pq* > 0

Therefore, one gets that 1 + *a*<sup>2</sup> > *a*1.

2. Since *a*<sup>1</sup> > 0, it is sufficient to prove that 1 + *a*<sup>2</sup> > 0. By hypothesis 1 > *R*0, then

$$\mu\_0 + q + q(h)\delta\_l\mu\_0 + q(h)^2\delta\_l(\mu\_0\mathbb{C} + \beta\Lambda) > q(h)^2\delta\_l(\mu\_0\mathbb{C} + \beta\Lambda) > q(h)^2\delta\_l N\beta\Lambda e^{-\delta\_{\tilde{l}}\Lambda}.$$
 Accordingly

$$\begin{split} \mu\_{0} + q + \varrho(h)\delta\_{I}(\mu\_{0} + \mu\_{0}\varrho(h)\mathbb{C} + \varrho(h)\beta\Lambda) &> \varrho(h)^{2}\delta\_{I}N\beta\Lambda e^{-\delta\_{I\_{\rm E}}\Delta} \Longleftrightarrow \\ \mu\_{0} + q + \varrho(h)\delta\_{I}q &> \varrho(h)^{2}\delta\_{I}N\beta\Lambda e^{-\delta\_{I\_{\rm E}}\Delta} \Longleftrightarrow \\ \mu\_{0} + (1 + \varrho(h)\delta\_{I})q &> \varrho(h)^{2}\delta\_{I}N\beta\Lambda e^{-\delta\_{I\_{\rm E}}\Delta} \Longleftrightarrow \\ 1 + \frac{\mu\_{0} - \varrho(h)^{2}\delta\_{I}N\beta\Lambda e^{-\delta\_{I\_{\rm E}}\Delta}}{pq} &> 0 \end{split}$$

3. Given that

$$
\mu\_0 - \varrho(h)^2 \delta\_I N \beta \Lambda e^{-\delta\_{I\_\mathbb{E}} \Lambda} < \mu\_0 + \mu\_0 \varrho(h) \mathbb{C} + \varrho(h) \beta \Lambda + \varrho(h) \delta\_I \eta\_\prime
$$

then

$$\frac{\mu\_0 - \varrho(h)^2 \beta e^{-\delta\_{l\_\mathcal{E}} \Lambda} \Lambda N \delta\_I}{pq} < 1,$$

that is *a*<sup>2</sup> < 1.

Next, by virtue of Lemma 1 we have that the eigenvalues of the Jacobian matrix of the system (35) evaluated at *Pf* satisfy |*λi*| < 1. for *i* = 1, 2. Then *Pf* is locally asymptotically stable if R<sup>0</sup> < 1.

#### *4.3. Global Stability of the NSFD Numerical Scheme*

Several authors have used discrete Lyapunov functions to study the global behavior of numerical solutions generated by non-standard finite difference schemes (*NSFD*), see [19,106–110]. For the study of the global stability of the critical points of the numerical scheme (34) it is necessary to use the Lyapunov functions and Equation (28).

**Theorem 8.** *The disease-free equilibrium point Pf* = 0, 0, <sup>Λ</sup> *μ*0 *of the scheme (34) is globally asymptotically stable if* R<sup>0</sup> ≤ 1*.*

**Proof.** Let the following Lyapunov function be

$$\mathcal{L}^{n} = \frac{T^{0}g\left(\frac{T^{n}}{T^{0}}\right) + \epsilon^{\delta\_{I\_{E}}\Delta}I^{n} + \beta T^{0}\left(\frac{1}{\mathcal{C}} + \varphi(h)\right)V^{n}}{\varphi(h)} + \sum\_{i=n-m\_{1}}^{n-1} \beta T^{i+1}V^{i}.\tag{37}$$

Using the inequality ln *<sup>z</sup>* <sup>≤</sup> *<sup>z</sup>* <sup>−</sup> 1, the difference of <sup>L</sup>*<sup>n</sup>* in (37) satisfies

$$\begin{split} \mathcal{L}^{n+1} - \mathcal{L}^{n} &\leq \left( \frac{T^{n+1} - T^{0}}{T^{n+1}} \right) \Big( \Lambda - \beta T^{n+1} V^{n} - \mu\_{0} T^{n+1} \Big) \\ &+ \epsilon^{\delta\_{l\_{\mathbb{E}}} \Delta} \Big( \beta T^{n-m\_{1}+1} V^{n-m\_{1}} e^{-\delta\_{l\_{\mathbb{E}}} \Delta} - \delta\_{l} I^{n+1} \Big) \\ &+ \frac{\beta T^{0}}{\mathcal{C}} \Big( N \delta\_{I} I^{n+1} - \mathcal{C} V^{n+1} - \beta T^{n} V^{n+1} \Big) + \beta T^{0} \Big( V^{n+1} - V^{n} \Big) \\ &+ \beta T^{n+1} V^{n} - \beta T^{n-m\_{1}+1} V^{n-m\_{1}} \\ &\leq -\mu\_{0} \frac{T^{n+1} - T^{0}}{T^{n+1}} - \epsilon^{\delta\_{l\_{\mathbb{E}}} \Delta} \delta\_{l} I^{n+1} (1 - \mathcal{R}\_{0}) . \end{split}$$

It follows that <sup>L</sup>*n*+<sup>1</sup> − L*<sup>n</sup>* <sup>≤</sup> 0 for all *<sup>n</sup>* <sup>≤</sup> 0 if <sup>R</sup><sup>0</sup> <sup>≤</sup> 1. This means that <sup>L</sup>*<sup>n</sup>* is monotone decreasing sequence and since <sup>L</sup>*<sup>n</sup>* <sup>≥</sup> 0 for *<sup>n</sup>* <sup>≥</sup> 0 then there exists a limit, i.e., lim*n*→<sup>∞</sup> <sup>L</sup>*<sup>n</sup>* <sup>≥</sup> 0. Hence lim*n*→<sup>∞</sup> <sup>L</sup>*n*+<sup>1</sup> − L*n* <sup>=</sup> 0, which implies that lim*n*→<sup>∞</sup> *<sup>T</sup><sup>n</sup>* <sup>=</sup> *<sup>T</sup>*0, lim*n*→<sup>∞</sup> *<sup>I</sup><sup>n</sup>* <sup>=</sup> 0, and from (34) lim*n*→<sup>∞</sup> *<sup>V</sup><sup>n</sup>* <sup>=</sup> 0. By the previous analysis, we conclude that *Pf* is globally asymptotically stable for scheme (34).

## **5. Numerical Simulations Using the** *NSFD* **Scheme**

In this section, we present some numerical solutions of the mathematical model of HIV. To carry out the simulations we use the constructed *NSFD* numerical scheme (34). We choose the parameter values based on existing experimental data and previous model studies, see [54,106,111,112]. The values of these parameters are given in Tables 1 and 2. For the first case we choose values such R<sup>0</sup> < 1 and for the second we have the case R<sup>0</sup> > 1.

**Table 1.** Parameter values for the numerical simulations when R<sup>0</sup> < 1.



**Table 2.** Parameter values for the numerical simulations when R<sup>0</sup> > 1.

For the numerical simulations we use a small step size, *h* = 0.001. For the discrete derivatives given in system (33), we have many options for the denominator function *<sup>ϕ</sup>*. We have chosen *<sup>ϕ</sup>*(*h*)=(<sup>1</sup> <sup>−</sup> *exp*(−*hp*))/*p*, where *<sup>p</sup>* <sup>=</sup> max *A*, Δ, *μ*0, *δI*, *δIE* are parameters of the model and included in the numerical scheme (33). This particular *ϕ* usually provides better numerical results based on previous articles related to *NSFD* schemes [113,114]. In addition, this option satisfies the asymptotic relation *ϕ*(*h*) = *h* + *O*(*h*2), and Rule 1. We performed numerical simulations to show that the solutions obtained by the proposed *NSFD* scheme and the well-known MATLAB routine dde23 agree. A great advantage of the proposed *NSFD* numerical scheme (34) is that the computation time is smaller. For the first case, the numerical solution using the *NSFD* scheme is obtained in approximately 0.083992 s, while the dde23 routine spent 2.573003 s. For the second case, the numerical solution using the *NSFD* scheme is obtained in approximately 0.176461 s, while the routine dde23 spent 11.181812 s. The results obtained are shown in Figures 2 and 3 respectively. It can be seen that the numerical solution of the proposed *NSFD* numerical scheme (34) and the one obtained by means of the routine dde23 agree.

**Figure 2.** Simulation *NSFD* versus dde23 when R<sup>0</sup> < 1.

**Figure 3.** Simulation *NSFD* versus dde23 when R<sup>0</sup> > 1.

#### **6. Conclusions**

We proposed a mathematical model based on a set of delay differential equations that describe intracellular HIV infection. The model considers three different subpopulations of cells and the HIV virus. The mathematical model is formulated in such a way that takes into account the time between viral entry into a target cell and the production of new virions. Moreover, this time is included using a discrete time delay. We analyzed the local stability of the infection-free and endemic equilibrium states. By using a suitable Lyapunov functional and the LaSalle invariant principle, it is proved that if the basic reproduction ratio is less than unity, the infection-free equilibrium is globally asymptotically stable. In addition, we designed a non-standard difference scheme that preserves some properties of the continuous model. We prove that the constructed *NSFD* scheme has the same equilibrium points of the continuous model, and the disease-free equilibrium holds the same stability properties. As required by the constraints of the real phenomena, the solutions given by the numerical scheme satisfy positivity and boundedness. The numerical simulations corroborate that the developed *NSFD* numerical scheme preserves the properties of the continuous model and presents a robust behavior when working with a variety of parameter values.

**Author Contributions:** Conceptualization, G.G.-P.; Formal analysis, G.G.-P., A.J.A. and J.J.N.; Investigation, G.G.-P. and A.J.A.; Methodology, G.G.-P., A.J.A., M.C. and N.D.L.E.; Software, G.G.-P.; Supervision, G.G.-P. and A.J.A.; Validation, A.J.A.; Visualization, A.J.A., J.J.N., M.C. and N.D.L.E.; Writing—original draft, G.G.-P.; Writing—review and editing, G.G.-P., M.C. and N.D.L.E. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors are grateful to the reviewers for their careful reading of this manuscript and their useful comments to improve the content of this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

