We discuss a quartic eigenvalue problem arising in the context of an optical waveguiding problem involving atomically thick 2D materials. The waveguide configuration we consider consists of a gradient-index (spatially dependent) dielectric equipped with conducting interior interfaces... This leads to a quartic eigenvalue problem with mixed transverse electric and transverse magnetic modes, and strongly coupled electric and magnetic fields. We derive a weak formulation of the quartic eigenvalue problem and introduce a numerical solver based on a quadratification approach in which the quartic eigenvalue problem is transformed to a spectrally equivalent companion problem. We verify our numerical framework against analytical solutions for prototypical geometries. As a practical example, we demonstrate how an improved quality factor (defined by the ratio of the real and the imaginary part of the computed eigenvalues) can be obtained for a family of gradient-index host materials with internal conducting interfaces. We outline how this result lays the groundwork for solving related shape optimization problems. (read more)