For NWs obeying a charge-neutrality constraint, Eq. ( ) is used to calculate ϵ F . The find_epsilon_F.m m-file uses the last energy level previously calculated by pdeeig as an initial guess to the built-in MATLAB root-finding fzero function. The fzero command solves the charge neutrality condition in charge_neutral.m, which calculates the difference between the integrated electron density (Eq. ( )) and the integrated depletion region (cf. Fig. ). If the ϵ F value predicted by find_epsilon_F.m is higher than the last energy level previously calculated by pdeeig, the lower and upper limit of the r vector is incrementally changed and inputted to pdeeig until all the required wavefunctions/energies are obtained and Eq. ( ) is satisfied.