h= 6.626e-34;

c=3e8;

k=1.38e-23;

T=300;

syms v

rho=8*pi*h*v.^3/c^3./(exp(h*v/k/T)-1); %函式

drho=diff(rho,v); %導數

vm=double(solve(drho)); %求導數為0的點，也就是極值點

lm=c/vm; %波長值

%

x=linspace(1,1e14,500); %自變數x取值，也就是v

y=subs(rho,v,x); %計算函式值，也就是rho

rm=subs(rho,v,vm); %極值點函式值

plot(x,y);hold on; %畫函式影象

plot(vm,rm,"r*");hold off; %標出極值點位置

xlabel("\it\nu");ylabel("\it\rho");

text(vm,rm+1e-20,...

["\it\nu_m\rm=",num2str(vm,"%.2e"),...

" \it\lambda_m\rm=",num2str(lm,"%.2e"),...

" \it\rho_m\rm=",num2str(rm,"%.2e")]);