I am not certain what the problem is.
The approach you took appears to work here (R2023b) —
c = 3e+08; % Speed of light in m/s
h = 6.63e-34; % Planck's constant in Js
k = 1.38e-23; % Boltzmann's constant in J/K
T = 5973.15; % Temperature in K
y = 0:10e-9:2000e-9; % Wavelength in m
R = ((2*pi*h*c^2)./y.^5).*(1./(exp((h*c)./(k*T*y))-1)); % Spectral Irridiance in W/m^2/m
plot(y,R)
xlabel({'Wavelength','(nm)'});
ylabel({'Spectral Irridiance','(W/m^2/nm)'});
title('Black Body Radiation Spectrum at 5700 Degrees Celcius');
ax = gca;
ax.XAxis.Exponent = -9;
.
