“…It is well known that exact analytical solutions to such problems are only available for simple structures such as a square or parabolic well [28], and even in these structures, in general, in the presence of perturbations such as external fields, disorder effects [29], etc., the problem cannot be solved exactly. There have been various numerical methods used to calculate the band profiles in QWs: the matrix approach (MA) [30,31], the transfer matrix (TM) method [32,33], the finite difference method (FDM) [34,35], the "shooting method" (SM) [34,36], the finite element (FE) technique [37,38], discreet variable representation (DVR) approach [39], envelope function (EF) method [40], Wentzel-Kramers-Brillouin (WKB) approximation [41], variational method (VM) [42], and Monte Carlo (MC) simulations [43]. Among them, the WKB and EF methods adopt approximations, thus giving unreliable results; the VM only works well at simple QWs and weak fields; the MC and FE methods are highly computerorientated approaches; the MA usually require wave function to be well behaved; the SM's speed comes at the cost of stability.…”