A formulation for the theoretical and numerical modeling of electromagnetic wave propagation in graphene-comprising waveguides is presented, targeting applications in the linear and nonlinear regime. Waveguide eigenmodes are rigorously calculated using the finite-element method (FEM) in the linear regime and are subsequently used to extract nonlinear properties in terms of the nonlinear Schrödinger equation framework. Graphene sheets are naturally represented as sheet/2D media and are seamlessly implemented with interface conditions in the FEM, thus greatly enhancing the computational efficiency. This formulation is used to analyze the nonlinear performance of several graphene-comprising waveguide configurations in the optical band, including silicon-based photonic waveguides, metal-based plasmonic waveguides and glass microfibers. Optimal design choices are identified for each configuration and subtle aspects of the FEM-based modeling, especially important for plasmonic waveguides, are highlighted.