*4.1. Equilibrium Points of the NSFD Numerical Scheme*

The equilibrium points of the scheme (34) are given by analyzing the behavior of system when *n* approaches to infinity. Thus, after a few calculations we find that

$$I^\* = \frac{\varrho(h)\beta T^\* V^\* e^{-\delta\_{I\_E}\Lambda} + I^\*}{1 + \varrho(h)\delta\_I} \tag{35}$$

$$\begin{split} V^\* &= \frac{\varrho(h)N\delta\_I I^\* + V^\*}{1 + \varrho(h)(C + \beta T^\*)} \\ T^\* &= \frac{\varrho(h)\Lambda + T^\*}{1 + \varrho(h)(\beta V^\* + \mu\_0)}. \end{split}$$

Note that the equations of the scheme (35) correspond to Equations (14)–(16). Thus, the critical points of the discrete scheme will coincide in the limit *h* → 0, with those of the continuous model.

#### *4.2. Local Stability of the NSFD Numerical Scheme*

For the study of the local stability of the critical points of the numerical scheme (34) it is necessary to use the following lemma:

**Lemma 1.** *The roots of the quadratic polynomial <sup>λ</sup>*<sup>2</sup> <sup>−</sup> *<sup>a</sup>*1*<sup>λ</sup>* <sup>+</sup> *<sup>a</sup>*<sup>2</sup> <sup>=</sup> 0, *satisfy* <sup>|</sup>*λi*<sup>|</sup> <sup>&</sup>lt; <sup>1</sup> *for <sup>i</sup>* <sup>=</sup> 1, 2 *if and only if the following conditions hold:*

*i.* 1 − *a*<sup>1</sup> + *a*<sup>2</sup> > 0, *ii.* 1 + *a*<sup>1</sup> + *a*<sup>2</sup> > 0, *iii. a*<sup>2</sup> < 1.
