Issues with self-consistent Poisson-Schrodinger solver

Question

I'm currently in the process of writing a self-consistent Schrodinger-Poisson solver for a device heterstructure (High Mobility Electron Transistor). The algorithm is based off of this journal⁽¹⁾. I am having a few difficulties however. It seems like I can either have Dirichlet BCs or Neumann BCs at the substrate (i.e. the right end of the structure). Below is a snapshot of a Dirichlet and Neumann BC, respectively where $ x_n=0 $ for the first, and $\partial x/\partial n = 0$ for the second plot. Note that the boundary in question $\partial \Omega$ is at the right end of the plot.

I'm really interested in having the boundary $\partial\Omega $ be like the plot below.

_{(1) A self‐consistent solution of Schrödinger–Poisson equations using a nonuniform mesh. Tan, I‐H. and Snider, G. L. and Chang, L. D. and Hu, E. L., Journal of Applied Physics, 68, 4071-4076 (1990)}

---

Below outlines the discretization scheme that I am using.

1. My plan is to simulate an AlGaAs/GaAs HEMT stack that looks almost identical to the one in the journal that I referenced above (The stack is in Fig. 3). There's only a few Angstroms of difference in some of my layers. I begin with a simple trial potential $$ V_{trial}(x) = -q\phi(x) + \Delta E_c(x) $$ Initially, this is a simple line that begins at 0.7 eV with a step function of 0.2 eV that spans the AlGaAs region.
2. I then solve the Schrodinger equation using the finite difference scheme.

The function Schrodinger($V_{trial}$) gets called first.

function [prob_Density, eig_Vals] = schrodinger(phi)
p = parameters;
n = p.grid_points;
t0 = p.hbar^2 / (2*p.m_AlGaN*p.dz^2);

A = zeros(n,n);

% Define first and last rows
A(1,1) = 2*t0; A(1,2) = -2*t0; A(n-1,n) = -2*t0; A(n,n) = 2*t0;      

for i=2:(n-1)  % Create Hamiltonian
    A(i,i) = 2*t0 + phi(i);
    A(i+1,i) = -t0; A(i,i+1) = -t0; 
    A(i-1,i) = -t0; A(i,i-1) = -t0;
end

[eig_states, eig_Vals] = eig(A);  % Solve for eigensolutions

prob_Density = eig_states.*conj(eig_states);  % Calculate probability function

end  % End of function

3. Using the result "prob_Density" from step 2, I then solve for the current distributions using $$n(x) = \sum_{k=1}^{m} \psi_k^* (x) \psi_k (x) n_k$$ where $$n_k = \frac{m^*}{\pi \hbar^2} \int_{E_k}^{\infty} \frac{1}{1+e^{(E-E_F)/kT}} dE $$

Next, I use the $n(x)$ and solve Poisson's equation. I once again use the finite difference scheme with the Poisson's equation. My Dirichlet boundary condition for x=0 is defined at the Schottky barrier height. My second boundary is of the Neumann type, where the derivative is zero.

The code for the "Charge Density" (i.e. $n(x)$) and Poisson solver are written below, respectively.

%#### FIND CHARGE DISTRIBUTION ########
function charge_density = n_x(prob_Density, eig_Vals, Ef, final_Eigstate)

p = parameters;

n = p.grid_points;
n_array = zeros(n,1);
const = (p.m_AlGaAs ./ (pi*p.hbar^2));

fermi_dirac = @(E) 1./(1+exp((E-Ef)./p.Vt)); % Unit eV

for eig_state = 1:final_Eigstate

    Ek = eig_Vals(eig_state,eig_state); % First,second,third .. eigenenergy
    nk = const .* integral(fermi_dirac, Ek, inf)*p.q;   %  [m]^-2
    n_x_arr = prob_Density(:,eig_state).*nk;    % [m]^-1 * [m]^-2 -> [m]^-3
    n_array = n_array + n_x_arr;
end

charge_density = n_array; % Electron volumetric density [m]^-3


end % End program


%######### START POISSON SOLVER ###########
function phi = poisson(n_x)

first_boundary = 0.7; % eV
second_boundary = 0;    

p=parameters;
n = p.grid_points;

% Create Poisson matrix elements 
A = zeros(n,n);

for i = 2:n-1
   A(i,i) = -2; 
   A(i+1,i) = 1; 
   A(i-1,i) = 1; 
   A(i,i+1) = 1; 
   A(i,i-1) = 1;
end

A(1,1) = -2; 
A(1,2) = 1;
A(n,n) = -1; % -1 a result of Neumann BC
A(n,n-1) = 1;

F = zeros(n,1);

%% RHS of matrix equation (the 'B' n x 1 term) AU=B, U=A\B
h = p.dz;

for i = 2:n-1 
   f = p.q*(p.ND(i) - n_x(i))./p.eps_hemt(i);
   F(i) = f*h^2;
end

f1 = p.q*(p.ND(1) - n_x(1))./p.eps_hemt(1);
fn = p.q*(p.ND(n) - n_x(n))./p.eps_hemt(n);

F(1) = f1*h^2 - first_boundary;

F(n) = fn*h^2 - second_boundary;


phi = (A\F);                      %  Voltage profile

end
%##### END POISSON SOLVER

4. Finally, I use the result from the Poisson solver and return to Step 2. I once again run the function Schrodinger($V_{new}$) and iterate through this process until a self consistent solution is found. This method brings me to the very first plot you see above, and I'm aiming to simulate the third plot.