A numerical method for calculating cutoff frequencies of modes of weakly guiding optical fibers is described. Starting from a properly constructed field representation, the method follows a variational approach to obtain an eigenvalue problem for the cutoff frequencies which can be easily solved by a standard numerical routine. The field representation uses an expansion in terms of the Laguerre-Gauss functions in the fiber core and limiting form for small arguments of modified Bessel functions in the fiber cladding. Fibers with power-law refractive index profiles and fibers with a profile showing a dip on the axis are analyzed. Results obtained for modes of arbitrary order (both azimuthal and radial) show that our method is very efficient and accurate.