黑体（将所有落到它上面的辐射全都吸收的系统）发出黑体辐射，太阳的辐射也差不多是黑体辐射。

clear;clc;
h=6.6262*10^(-34);%普朗克常数
c=2.997925*10^8;%光速
k=1.3807*10^(-23);%玻尔兹曼常数
T=5700;%太阳的温度
L=linspace(1e-8,1e-2,10000000);%设置波长范围
f1=8*pi*k*T./L.^4;
semilogx(L,f1)
ylim([0,1.2*1e+6])
title('黑体辐射的能量密度与波长的关系')
xlabel('波长（m=10^6um）')
ylabel('能量密度')
hold on
f2=8*pi*h*c*L.^(-5)./(exp(1).^(h*c./(L.*k*T))-1);
semilogx(L,f2)
xlim([0,1e+6])
legend('瑞利-吉恩斯定律','普朗克定律')

clear;clc;
h=6.6262*10^(-34);%普朗克常数
c=2.997925*10^8;%光速
k=1.3807*10^(-23);%玻尔兹曼常数
L=linspace(1e-8,1e-2,10000000);%设置波长范围
subplot(2,2,1)
T1=3400;
f1=8*pi*h*c*L.^(-5)./(exp(1).^(h*c./(L.*k*T1))-1);
semilogx(L,f1)
ylim([0,1e+6])
title('一等星（T=3400K）')
subplot(2,2,2)
T2=6400;
f2=8*pi*h*c*L.^(-5)./(exp(1).^(h*c./(L.*k*T2))-1);
semilogx(L,f2)
title('南河三（T=6400K）')
subplot(2,2,3)
T3=9200;
f3=8*pi*h*c*L.^(-5)./(exp(1).^(h*c./(L.*k*T3))-1);
semilogx(L,f3)
title('天狼星（T=9200K）')
subplot(2,2,4)
T=5700;
f=8*pi*h*c*L.^(-5)./(exp(1).^(h*c./(L.*k*T))-1);
semilogx(L,f)
title('太阳（T=5700K）')

用liveeditor有数值滑块，多移T就能出现一系列的曲线了。缺点就是太慢。
意外收获
mathpix截图就能获取LaTex公式，自己再稍微修改下，然后在matlab中插入即可。