Summary: Interface problems have long been a major focus of scientific computing, leading to the development of various numerical methods. Traditional mesh-based methods often employ time-consuming body-fitted meshes with standard discretization schemes or unfitted meshes with tailored schemes to achieve controllable accuracy and convergence rate. Along another line, mesh-free methods bypass mesh generation but lack robustness in terms of convergence and accuracy due to the low regularity of solutions. In this study, we propose a novel method for solving interface problems within the framework of the random feature method (RFM). This approach utilizes random feature functions in conjunction with a partition of unity as approximation functions, and solves a linear least-squares system to obtain the approximate solution. In the context of interface problems, two innovative and crucial components are incorporated into the RFM. Firstly, we utilize two sets of random feature functions on each side of the interface, allowing for the inclusion of low regularity or even discontinuous behaviors in the solution. Secondly, the construction of the loss function is based on the assessment of the partial differential equation, initial/boundary conditions, and the interface condition on collocation points. This approach ensures that these conditions are equally satisfied. Consequently, the challenges arising from geometric complexity primarily manifest in the generation of collocation points, a task amenable to standard methods. Importantly, the proposed method retains its meshfree characteristics and robustness when addressing problems featuring intricate geometries. We validate our method through a series of linear interface problems with increasingly complex geometries, including two-dimensional elliptic and three-dimensional Stokes interface problems, a three-dimensional elasticity interface problem, a moving interface problem with topological change, a dynamic interface problem with large deformation, and a linear fluid-solid interaction problem with complex geometry. Our findings show that despite the solution often being only continuous or even discontinuous, our method not only eliminates the need for mesh generation but also maintains high accuracy, akin to the spectral collocation method for smooth solutions. Remarkably, for the same accuracy requirement, our method requires two to three orders of magnitude fewer degrees of freedom than traditional methods, demonstrating its significant potential for solving interface problems with complex geometries and predetermined intricate evolution.