It’s well known that the barycentric form of the Lagrange interpolation polynomial is numerically stable and used to reconstruct complex functions from discrete data. The interpolation rational function \(\mathcal{R}_{m,n}\) such that \( \mathcal{R}_{m,n}(z_j)=\frac{\mathcal{P}_m(z_j)}{\mathcal{Q}_n(z_j)}=f_j\) for \(j=0,\dots,N\), with \(f_j=f(z_j)\) and \(z_i \neq z_j\), \(i\neq j\) known, and with nonnegative \(m\) and \(n\) satisfying \(N=m+n\) is the second approach in the interpolation problem. These methods become inefficient for functions with a branch cut or essential singularity. Suppose that \(a_j \in \operatorname{int} D_R\), \(j = 1, 2, \dots, m\), \(D_R=\{z: |z|<R\},\) are the singularities of \(f(z)\). The new barycentric rational interpolation function of the form \[ \mathcal{R}_{\mathbf{N}}(z)=\frac{\sum_{j=0}^m \frac{\sigma_j}{N_j}\sum_{k=0}^{N_j-1}\frac{z_{jk}-a_j}{z_{jk}-a}f_{jk}}{\sum_{j=0}^m\frac{\sigma_j}{N_j}\sum_{k=0}^{N_j-1}\frac{z_{jk}-a_j}{z_{jk}-a}} \] is proposed in this paper. This interpolation function is pole-free, exponentially convergent, and numerically stable, requiring only \(\mathcal{O}(N)\) operations. A convergence analysis of the method is provided. Numerical examples that illustrate the theoretical results and demonstrate the accuracy and efficiency of the methodology are given. Some applications of the method, including the numerical solutions of boundary value problems and the zero locations of a holomorphic functions, are discussed in the paper. The article contains many drawings, graphs and Matlab codes illustrating the considered methods.