With the wavefunctions and Fermi level computed, the summation n |ψ n ( x, ỹ)| 2 √ ϵ F -ϵ n is computed by the psi_sqrt_eps_summation.m m-file and stored in the global variable psi_sqrt_eps_sum (which is used by the n_e_func.m m-file to compute S e in Eq. ( )). It is worth mentioning that n_e_func.m utilizes the griddata_NaN.m and symmetrize_coordinates.m m-files to symmetrize the spatial electron density for the triangular and hexagonal geometries based on their respective two-and six-fold symmetries. The depletion region (cf. Fig. ) is computed and stored in the global variable heaviside_n_D, which contains a vector having elements of 1 for ( x, ỹ) values where V T ( x, ỹ) ≥ E F and 0 otherwise. For hexagonal geometries that do not contain polarization effects, the S D ( x, ỹ) and S e ( x, ỹ) source terms are computed in the n_D_func.m and n_e_func.m m-files, respectively, and their sum is used as input for the "f" coefficient in assempde to solve the dimensionless Poisson equation. For both of the Ga-and N-face triangular geometries, the sum of n_D_func.m, n_e_func.m, as well as the polarization source terms computed in rho_semipolar_inside.m and rho_strain_analytic.m are used as input for the "f" coefficient in assempde. The inputs for the "c" coefficient and boundary condition matrix, b, in assempde were discussed in Section 3.2.