Many physical systems can be described in terms of matrix models that we often cannot solve analytically. Fortunately, they can be studied numerically in a straightforward way. Many commonly used algorithms follow the Monte Carlo method, which is efficient for small matrix sizes but cannot guarantee ergodicity when working with large ones. In this paper, we propose an improvement of the algorithm that, for a large class of matrix models, allows to tunnel between various vacua in a proficient way, where sign change of eigenvalues is proposed externally. We test the method on two models: the pure potential matrix model and the scalar field theory on the fuzzy sphere.