Here's some code to get you started
% Data
T=300;
k=1.38064852e-23;
Ea=7.2*10^(-21);
Ed=7.2*10^(-21);
Nv=2.5*10^(25)*((0.59*T/300)^(3/2))*10^(-6);
Nc=2.5*10^(25)*((1.08*T/300)^(3/2))*10^(-6);
Ec=0.5*(1.166-0.000473*T*T/(636+T))*1.6*10^(-19);
Ev=0;
% Collect data to be passed to function
data = [T, k , Ea, Ed, Nv, Nc, Ec, Ev];
% Dopant concentrations
Nd = 10^17;
Na = 10^5;
Ef0 = 10^-21; % Initial guess at fermi energy
% Use fzero to find fermi energy, i.e. the value of Ef that makes
% function Efn return zero
Ef = fzero(@Efn, Ef0,[],data,Nd,Na);
disp(Ef)
function F = Efn(Ef,data, Nd, Na)
T = data(1);
k = data(2);
Ea = data(3);
Ed = data(4);
Nv = data(5);
Nc = data(6);
Ec = data(7);
Ev = data(8);
kT = k*T;
F = Nc*exp(-(Ec-Ef)/kT) + Na/(1+4*exp(-(Ef-Ea)/kT)) - Nv*exp(-(Ef-Ev)/kT) - Nd/(1+2*exp(-(Ed-Ef)/kT));
end
The above will calculate the Fermi level for one pair of dopant concentrations. See if you can take it from here.
