clear all;clc
format long

p=0.01:1:100;nni=4e15*p;
% 超越方程 nf/ni=1+6.94e25/ni*log(nf/ni)
for i=1:length(nni);
     ni=nni(i);
     [nf(i) fval(i)]=fsolve(@(x)x/ni-1-6.94e25/ni*log(x/ni),ni);
end;
[4e15*p ; nf ]'
fval'

figure(1)
plot(4e15*p,nf)