The critical challenge facing the modeling accuracy of borehole-to-surface electrical potential (BSEP) is the singularity effect resulted from the field source and the truncation error due to improper boundary conditions. However, different boundary conditions can reduce the truncation error by enlarging the computational domain. This enlargement undoubtedly causes a significant increase in the calculation burden. To some extent, the mixed boundary condition-based finite element method (MBC-FEM) can mitigate this problem. However, it is inapplicable to the BSEP modeling with line source because these algorithms were initially derived from the point source. For improving the BSEP modeling efficiency and precisely represent the specific correlation between BSEP observation results and the spatial location of the abnormal underground body, this paper introduces a novel finite-infinite coupling method (FICM) based on the secondary potential of arbitrary line source for three-dimensional BSEP modeling. By employing the variational principle and the different interpolation shape functions of both the FEM and the infinite element method (IEM), we discretize the function of three-dimensional line field source in the elements. The proposed coupled algorithm is then tested on a uniform half-space model, a low-resistance spherical anomaly model and an international standard model. The results show the proposed FICM can eliminate the singularity effect to reduce the truncation error. More importantly, the calculation efficiency is significantly improved by essential virtue of the needlessness of enlarging the computational domain, as well as the one-off generation of coefficient matrices. In this sense, the proposed FICM represents a more appropriate candidate for BSEP forward and inversion when a large-scale calculation is involved, especially in terms of irregular underground abnormal bodies in the line source conditions.