ND = 1e16;

for ii = 1:20;
    if ii == 1
        
        NDp (1) = ND 
    else      
        % Calculate Eg
        T=77;
        Eg = 1.169-4.9e-4*T^2/(T+655);

        % calculate KT

        KT= 0.0259*T/300;

        %for Phosphorous
        Ed = 0.046;

        %Calculate ND+
      
        ED = Eg - Ed; 

        % Ef we don't know. Let's iterate
        NC = 2.9e19*((77/300)^(3/2));

        EF (ii) = Eg + KT*log(NDp(ii-1)/NC)

        NDp (ii) = ND/(1+2*exp((EF(ii)-ED)/KT))
    end
end

figure; plot(NDp)
figure; plot(EF)