  The proximity to the Eddington luminosity has been attributed as the cause of
several observed effects in massive stars. Computationally, if the luminosity
carried through radiation exceeds the local Eddington luminosity in the
low-density envelopes of massive stars, it can result in numerical
difficulties, inhibiting further computation of stellar models. This problem is
exacerbated by the fact that very few massive stars are observed beyond the
Humphreys-Davidson limit, the same region in the Hertzsprung-Russell diagram
where the aforementioned numerical issues relating to the Eddington luminosity
occur in stellar models. One-dimensional stellar evolution codes have to use
pragmatic solutions to evolve massive stars through this computationally
difficult phase. In this work, we quantify the impact of these solutions on the
evolutionary properties of massive stars. We used the stellar evolution code
MESA with commonly used input parameters for massive stellar models to compute
the evolution of stars in the initial mass range of 10-110 M$_\odot$ at
one-tenth of solar metallicity. We find that numerical difficulties in stellar
models with initial masses greater than or equal to 30 M$_\odot$ cause these
models to fail before the end of core helium burning. Recomputing these models
using the same physical inputs but three different pragmatic solutions to treat
the numerical instability, we find that the maximum radial expansion achieved
by stars can vary by up to 2000 R$_\odot$, while the remnant mass of the stars
can vary by up to 14 M$_\odot$ between the sets. These differences can have
implications on studies such as binary population synthesis.
