希望能够绘制出如图片所示的效果

以下是MATLAB的代码：

%温度不同的普朗克黑体单色辐射能力与波长的曲线
clear %清除变量
k=1.38054e-23; %玻尔兹曼常数
h=6.626e-34; %普朗克常数
c=2.997925e8; %光速
sigma=5.6688e-008; %斯特潘常数
b=0.0029; %维恩常数
t=1400:100:2000; %热力学温度向量
n=length(t); %向量长度
lambda=[0:0.01:5]*1e-6; %波长向量
lambda(1)=eps; %给零加一小量使分母不为零
[T,L]=meshgrid(t,lambda); %波长和温度矩阵
M=2*pi*h*c^2./(exp(h*c./(k*T.*L))-1)./L.^5;%单色辐射能力,14.1.2
figure %创建图形窗口
plot(lambda*1e6,M) %画曲线

hl=legend([repmat('\itT\rm=',n,1),num2str(t'),repmat('K',n,1)]);%标记图例
fs=16; %字体大小
set(hl,'fontsize',fs) %设置图例大小
grid on %加网格
title('普朗克黑体单色辐射能力与波长的关系','fontsize',fs)%标题
xlabel('波长\it\lambda\rm/\mum','fontsize',fs)%横坐标
yl='单色辐射能力\itM\rm(\it\lambda\rm,\itT\rm)/(W\cdotm^-^3)';%纵坐标字符串
ylabel(yl,'fontsize',fs) %纵坐标
txt=['\itb\rm=' num2str(b) 'm\cdotK']; %维恩常数文本
txt=[txt ',\it\sigma\rm=' num2str(sigma) 'W/(m^2\cdotK^4)'];%斯特潘常数文本
text(0,max(M(:))/10,txt,'fontsize',fs) %显示常数
hold on %保持图像
[mx,ix]=max(M); %找最大值和下标
stem(lambda(ix)*1e6,mx,'--','filled') %画直杆图
text(lambda(ix)*1e6,mx,[num2str(lambda(ix)'*1e6)],'fontsize',fs)%显示峰值波长
t=1300:2020; %较密的温度向量
lambda=b./t; %波长向量
m=2*pi*h*c^2./(exp(h*c./(k*t.*lambda))-1)./lambda.^5;%单色辐射能力向量
plot(lambda*1e6,m) %画峰值曲线

亲，您好，我是CSDN必问的Q妹，你这个问题在必问区提问可能会更快速解决哦，邀您体验：https://biwen.csdn.net/

这是我自己改写的，有些不完善，你可以再改一下。

encoding: utf-8

import matplotlib.pyplot as plt
from sympy import *

import numpy as np

k=1.38054e-23

h=6.626e-34

c=2.997925e8

sigma=5.6688e-008

b=0.0029
t1 = 1400
t2 = 1500
t3 = 1600
t4 = 1700
t5 = 1800
t6 = 1900
t7 = 2000
t8 = 6000
array

x=list((x*0.01+1e-10)*1e-6 for x in range(0,500))
m=np.zeros(500)

def f(a):
for i in range(0,500):
m[i]=2*np.pi*h*c**2/(exp(h*c/(k*a*x[i]))-1)/x[i]**5
return m

g = [x[np.argmax(f(t1))], x[np.argmax(f(t2))], x[np.argmax(f(t3))], x[np.argmax(f(t4))], x[np.argmax(f(t5))], x[np.argmax(f(t6))] , x[np.argmax(f(t7))]]
n = [max(f(t1)) , max(f(t2)) ,max(f(t3)) ,max(f(t4)) ,max(f(t5)) ,max(f(t6)) ,max(f(t7))]

print(np.argmax(f(t1)) , max(f(t1)))
print(np.argmax(f(t2)) , max(f(t2)))
print(np.argmax(f(t3)) , max(f(t3)))
print(np.argmax(f(t4)) , max(f(t4)))
print(np.argmax(f(t5)) , max(f(t5)))
print(np.argmax(f(t6)) , max(f(t6)))
print(np.argmax(f(t7)) , max(f(t7)))

plt.plot(x,f(t1),label='T=1400K')
plt.plot(x,f(t2),label='T=1500K')
plt.plot(x,f(t3),label='T=1600K')
plt.plot(x,f(t4),label='T=1700K')
plt.plot(x,f(t5),label='T=1800K')
plt.plot(x,f(t6),label='T=1900K')
plt.plot(x,f(t7),label='T=2000K')
plt.plot(g , n , color='red', linewidth=2, linestyle='--')
plt.stem(g , n , linefmt='b:' , markerfmt='C3.', basefmt='r-')

for a, b in zip(g, n):
plt.text(a, b, a, ha='left', va='bottom', fontsize=10)

plt.title(u'普朗克黑体单色辐射能力与波长的关系')
plt.xlabel(u'波长λ\m')
plt.ylabel(u'单色辐射能力M(λ,T)/(W▪m^)')

plt.legend()
plt.show()