For highly compact neutron stars , constructing numerical initial data for black hole–neutron star binary evolutions is very difficult . We describe improvements to an earlier method that enable it to handle these more challenging cases . We examine the case of a 6:1 mass ratio system in inspiral close to merger , where the star is governed by a polytropic \Gamma = 2 , an SLy , or an LS220 equation of state . In particular , we are able to obtain a solution with a realistic LS220 equation of state for a star with compactness \num 0.26 and mass \SI 1.98 M _ { \odot } , which is representative of the highest reliably determined neutron star masses . For the SLy equation of state , we can obtain solutions with a comparable compactness of 0.25 , while for a family of polytropic equations of state , we obtain solutions with compactness up to \num 0.21 , the largest compactness that is stable in this family . These compactness values are significantly higher than any previously published results . We find that improvements in adapting the computational domain to the neutron star surface and in accounting for the center of mass drift of the system are the key ingredients allowing us to obtain these solutions .