Population dynamics indicate the changes in size and composition of population through time, as well as biotic and abiotic factors influencing those changes. Predator-prey (PP) relationship with harvesting and functional response involving prey refuge with Holling type I functional response (HTIFR) has been studied with recommendations on their extension to include Holling type II functional response (HTIIFR). There persists a problem in fifinding the numerical solution of predator-prey system having HTIIFR particularly when optimal selective harvesting (OSH) is carried out. Therefore, this study focuses on numerical analysis of predator-prey model with HTIIFR having time delay OSH. The objectives of this study includes to: Formulate a HTIIFR PP model with time delay OSH; Carry outstability analysis of the model and; Simulate the numerical outcome of the model. The methodology involves a step by step formulation of the model by considering various parameters and useful assumptions. This is followed by carrying out stability analysis of the model through local and global stability analysis techniques. Lastly, MATLAB software isutilized using Runge-Kutta technique in simulating the numerical outcomes of the model. The results of the study includes a formulated model which is both LAS and GAS when γ1 ≥ 0.00604. The numerical analysis has shown that the rate harvesting goes down due to a rise in the cost of harvesting. This study is useful in informing policy makers on how to control interaction between predator-prey environments so that extinction is controlled.