The coupling of electrothermal and optical models in multidimensional semiconductor laser simulation is a non-trivial task: While electrothermal models are conveniently formulated in a system of partial differential equations, the optical problem necessitates the solution of Maxwell's equations with results in eigenvalue problems with the eigenvalues representing lasing wavelengths and losses of the laser's modes. The state-of-the-art approach to achieve a self-consistent electro-optothermal solution os a Gummel-type iteration where the electrothermal equations (using a Newton scheme), and the optical eigenvalue problems are solved iteratively until a convergence criterion is reached. This is extremely time-consuming and unstable, especially for simulation of devices featuring closely spaced multiple modes (for example broad area lasers, VCSELs). In this work, we present a novel numerical electro-optothermal coupling scheme for semiconductor lasers which is based on the integration of the optical eigenvalue into a single global Newton formulation, The necessary derivatives are obtained by perturbation theory. The proposed scheme is more robust and decreases the computational burden (simulation time) by more than an order of magnitude compared to Gummel-type interations. This novel coupling scheme allows to investigate the influence of the electro-optothermal coupling on the device characteristics and internal physics. The effects of bias conditions on the modal dynamics, optical wavelengths, losses, and far-field patterns are analyzed. For a VCSEL, we quantify the role of gain guiding and thermal lensing.