I'm working on a project that requires me to numerically integrate the Planck spectral distribution. The object is to find the median wavelength, with exactly 50% of the radiance on either side. I'm using a standard composite Simpson rule method and I get good convergence with temps around 5000K. For temps at and above 7000K, I can't get convergence. Is there a better scheme to use for this application?