The SIS model is one of the fundamental frameworks used to describe infectious diseases in which recovered individuals do not acquire permanent immunity. However, traditional numerical schemes such as the Euler and Runge–Kutta methods often fail to preserve key qualitative properties of the model—namely, the positivity of solutions, conservation of total population, and dynamic consistency—especially when large time steps are used. In this study, we propose a novel Nonstandard Finite Difference (NSFD) scheme for the SIS model. The scheme guarantees the positivity of solutions, conservation of total population, and elementary stability, while exhibiting weak dependence on both the step size and parameter variations. Theoretical analysis and numerical experiments demonstrate that the proposed scheme provides more stable, robust, and dynamically consistent solutions than the standard Euler, second-order Runge–Kutta (RK2), and Mickens-type NSFD methods. The results indicate that the proposed NSFD approach offers an accurate and reliable numerical tool for simulating SIS and other structurally similar epidemic models.