废话不多说，直接上图片
% lambda为0.85 % mu为0.85/0.3 % rho 为0.3 % R 为100 % h 为5 clc; clear; syms a T N k R b h a=0.85; R=100; h=5; b=0.3; N=4; f1=N*(N-1)*(1-symsum(exp(-a*T)*((a*T)^(k))/(factorial(k)),k,0,N-1)); f2=symsum(exp(-a*T)*((a*T)^(k))/(factorial(k-2)),k,2,N-1); f3=0; for m=1:N for i=0:m-1 f3=f3-exp(-a*T)*((a*T)^i)/factorial(i); end end f3=f3+N; f4=2*(exp(-a*T)+f3); f5=h*(f1+f2)/f4+R*(a*(1-b)/(exp(-a*T)+N-symsum(exp(-a*T)*(a*T)^(k)*(N-k)/(factorial(k)),k,0,N-1))); F=diff(f5,T); % T=vpasolve(F==0,T); f6=R*(a*(1-b)/(exp(-a*T)+N-symsum(exp(-a*T)*(a*T)^(k)*(N-k)/(factorial(k)),k,0,N-1))) F2=diff(f6,T) T=vpasolve(F2==0,T)