In this contribution, the complicated dynamical behaviors and optimal harvesting policy of a discrete-time predator–prey model with fear and refuge effects are formulated. Both the fear and prey refuge effects refer to an interaction between predator and prey. In the first place, the existence and local stability of three fixed points of proposed model are investigated by virtue of our methodology, that is, the eigenvalues of the Jacobian matrix. One step further, it is worth mentioning that the model undergoes flip bifurcation (i.e., period–doubling bifurcation) and Neimark–Sacker bifurcation at the interior fixed point by the utilization of bifurcation theory and center manifold theory. Also, optimal harvesting strategy is investigated, and the expressions of optimal harvesting efforts are determined. Two examples, in the end, are put forward to prove that they are consistent with the previous theoretical results.

In ecology, the most significant aspect is that the interactions between predators and prey are extremely complicated. Numerous experiments have shown that both direct predation and the fear induced in prey by the presence of predators lead to a reduction in prey density in predator-prey interactions. In addition, a suitable shelter can effectively stop predators from attacking as well as support the persistence of prey population. There has been less exploration of the effects of not only fear but also refuge factors on the dynamics of predator prey interactions. In this paper, we unveil several conclusions about a predator-prey system with fuzzy parameters, considering the cost of fear in two prey species and the effect of shelter on two prey species and one predator. As the first step of the investigation, the boundedness and non-negativity of the solutions to the system are put forward. Using the Jocabian matrix and Lyapunov function methods, we further analyze the existence and stability of the available equilibria and also the existence of Hopf bifurcation, considering the fear parameter as the bifurcation parameter that has been observed by applying the normal form theory. Finally, numerical simulations help us better understand the dynamics of the model, in which some interesting chaotic phenomena are also exhibited.