A general numerical algorithm for the computation of the fundamental modes and related cutoff wavenumbers in arbitrary shaped inhomogeneous waveguides is presented. The method exploits a generalized Donsker-Kac formula to express the lowest order modes in terms of asymptotic generalized Wiener-Ito integrals, whose computation is carried out by means of Monte Carlo methods. Comparison with known solutions and computational budget indicate that the proposed method is indeed accurate, versatile, as well as computationally efficient.