We present radiation hydrodynamics simulations of the collapse of massive pre-stellar cores . We treat frequency dependent radiative feedback from stellar evolution and accretion luminosity at a numerical resolution down to 1.27 AU . In the 2D approximation of axially symmetric simulations , it is possible for the first time to simulate the whole accretion phase ( up to the end of the accretion disk epoch ) for the forming massive star and to perform a broad scan of the parameter space . Our simulation series show evidently the necessity to incorporate the dust sublimation front to preserve the high shielding property of massive accretion disks . While confirming the upper mass limit of spherically symmetric accretion , our disk accretion models show a persistent high anisotropy of the corresponding thermal radiation field . This yields to the growth of the highest-mass stars ever formed in multi-dimensional radiation hydrodynamics simulations , far beyond the upper mass limit of spherical accretion . Non-axially symmetric effects are not necessary to sustain accretion . The radiation pressure launches a stable bipolar outflow , which grows in angle with time as presumed from observations . For an initial mass of the pre-stellar host core of 60 , 120 , 240 , and 480 \mbox { M } _ { \odot } the masses of the final stars formed in our simulations add up to 28.2 , 56.5 , 92.6 , and at least 137.2 \mbox { M } _ { \odot } respectively .