In this paper, we construct new nonstandard finite difference schemes to approximate a set of positive solutions for the predator–prey model, which contains different functional responses. The organization of the denominator of the discrete derivative and nonlocal approximations of nonlinear terms are employed to design the new schemes. The approach results in significant qualitative improvements in how the numerical solution behaves. We establish that the proposed nonstandard finite difference methods are elementary stable and satisfy the positivity requirement. In addition, the instances of applying PESN methods to some predator–prey systems using the Beddington–DeAngelis and Nicholson–Bailey functional responses are provided here. Finally, some numerical comparisons are presented to illustrate our findings. Our results indicate that the proposed methods are very suitable for the symmetric model of predator–prey.
