A new algorithm, denoted by RSRR, is presented for solving large-scale nonlinear eigenvalue problems (NEPs) with a focus on improving the robustness and reliability of the solution, which is a challenging task in computational science and engineering. The proposed algorithm utilizes the Rayleigh–Ritz procedure to compute all eigenvalues and the corresponding eigenvectors lying within a given contour in the complex plane. The main novelties are the following. First and foremost, the approximate eigenspace is constructed by using the values of the resolvent at a series of sampling points on the contour, which effectively circumvents the unreliability of previous schemes that using high-order contour moments of the resolvent. Secondly, an improved Sakurai–Sugiura algorithm is proposed to solve the projected NEPs with enhancements on reliability and accuracy. The user-defined probing matrix in the original algorithm is avoided and the number of eigenvalues is determined automatically by the provided strategies. Finally, by approximating the projected matrices with the Chebyshev interpolation technique, RSRR is further extended to solve NEPs in the boundary element method, which is typically difficult due to the densely populated matrices and high computational costs. The good performance of RSRR is demonstrated by a variety of benchmark examples and large-scale practical applications, with the degrees of freedom ranging from several hundred up to around one million. The algorithm is suitable for parallelization and easy to implement in conjunction with other programs and software.