In this paper, we study a multi-server queueing system with retrials and an infinite orbit. The arrival of primary customers is described by a batch Markovian arrival process (BMAP), and the service times have a phase-type (PH) distribution. Previously, in the literature, such a system was mainly considered under the strict assumption that the intervals between the repeated attempts from the orbit have an exponential distribution. Only a few publications dealt with retrial queueing systems with non-exponential inter-retrial times. These publications assumed either the rate of retrials is constant regardless of the number of customers in the orbit or this rate is constant when the number of orbital customers exceeds a certain threshold. Such assumptions essentially simplify the mathematical analysis of the system, but do not reflect the nature of the majority of real-life retrial processes. The main feature of the model under study is that we considered the classical retrial strategy under which the retrial rate is proportional to the number of orbital customers. However, in this case, the assumption of the non-exponential distribution of inter-retrial times leads to insurmountable computational difficulties. To overcome these difficulties, we supposed that inter-retrial times have a phase-type distribution if the number of customers in the orbit is less than or equal to some non-negative integer (threshold) and have an exponential distribution in the contrary case. By appropriately choosing the threshold, one can obtain a sufficiently accurate approximation of the system with a PH distribution of the inter-retrial times. Thus, the model under study takes into account the realistic nature of the retrial process and, at the same time, does not resort to restrictions such as a constant retrial rate or to rough truncation methods often applied to the analysis of retrial queueing systems with an infinite orbit. We describe the behavior of the system by a multi-dimensional Markov chain, derive the stability condition, and calculate the steady-state distribution and the main performance indicators of the system. We made sure numerically that there was a reasonable value of the threshold under which our model can be served as a good approximation of the BMAP/PH/N queueing system with the PH distribution of inter-retrial times. We also numerically compared the system under consideration with the corresponding queueing system having exponentially distributed inter-retrial times and saw that the latter is a poor approximation of the system with the PH distribution of inter-retrial times. We present a number of illustrative numerical examples to analyze the behavior of the system performance indicators depending on the system parameters, the variance of inter-retrial times, and the correlation in the input flow.