Setting things up

k = 1.381*^-23;
e = 1.602*^-19;
Eg  = 1.124 e; (* band gap of Si in eV, multiplying by e to convert to J *)
Ev=0;  (* picking my zero of energy to be at VB *)
Ec=Ev+Eg;
Ed = Ec-0.045 e;   (* the absolute position of donor level, not relative *)
Nd = 1*^21;

(* I'm going to multiply mu by charge of electron in these next equations so mu in the plot will have units of eV. I will also divide all concentrations by 10^21 for now so the "FindRoot" command doesn't have to work with huge exponents. *)
n[mu_,T_]:= 2.86*^25 (T/300)^(3/2) Exp[(mu e-Ec)/(k T)]/1*^21  
p[mu_,T_]:= 2.50434*^25 (T/300)^(3/2) Exp[(Ev - mu e)/(k T)]/1*^21  
Ndplus[mu_,T_]:=Nd/(1 + 2 Exp[(mu  e- Ed)/(k T)])/1*^21(* mu in eV, n in 10^21 *)

(* The extrinsic condition is "n = p + Ndplus". I'm going to plot both sides of that equation as a function of mu then find the mu where they intersect. Then I'll use that mu in the n equation to get the answer for a desired temperature. *)

Solving when T = 77K

T=77;
Plot[{n[mu,T], p[mu,T]+Ndplus[mu,T]},{mu,0, 1.2}]  (* plotting both sides of the extrinsic equation *)

(* Looks like answer must be around mu = 1.05 *)
FindRoot[n[mu,T]== p[mu,T]+Ndplus[mu,T],{mu,1.05}]

foundmu=mu /.%;  (* This is needed because the FindRoot command provides its answer in a weird format *)
n[foundmu,T]  * 1*^21 (* the answer for T=77! *)

Solving when T = 300K

(* I'll do everything in one cell now. I'll try to use the same starting point for my in the FindRoot command that I did for 77 K, namely 1.05. It could fail... but fortunately it doesn't. *)
T=300;
Plot[{n[mu,T], p[mu,T]+Ndplus[mu,T]},{mu,0, 1.2}]
FindRoot[n[mu,T]== p[mu,T]+Ndplus[mu,T],{mu,1.05}]  
foundmu=mu /.%;  
n[foundmu,T]  * 1*^21(* the answer for T=300! *)

Solving for arbitrary temperatures

Off[General::munfl]  (* This suppresses a warning message that will otherwise show up *)
(* I'm going to skip making the plot of the two sides of the extrinsic equation for the desired temperature and will skip straight to the FindRoot command *)
actualmu[T_]:=mu/. FindRoot[n[mu,T]== p[mu,T]+Ndplus[mu,T],{mu,1.05}];
Plot[actualmu[T],{T,5,1000},PlotRange->{0,1.2}]

actualn[T_]:=n[actualmu[T],T] * 10^21
LogPlot[actualn[T],{T,15,1000},PlotRange->{{0,1000},{10^19,10^23}}]