I am trying to solve the $\textbf{1-D}$ Poisson equation for a semiconductor structure at equilibrium (There is no external bias applied).
$\textbf{Background}$
\begin{equation} \frac{d^2V}{dx^2} = -\frac{\rho(V)}{\epsilon}\\ \rho(V) = q(N_D-N_A+p(V)-n(V)) \end{equation}
When there are two different semiconductors placed in contact, there will be a redistribution of charge, causing a dependence of charge on potential (the bands will bend).
$\textbf{Solving with the Newton-Rhapson method}$
I am trying to solve the above non-linear equation using the newton-rhapson method. Taking the central difference, the second derivative becomes: \begin{equation} \frac{d^2V}{dx^2}|_i = \frac{(V_{i+1}-2V_i+V_{i-1})}{dx^2} \end{equation}
The issue with solving this now is that this derivative cannot be computed for the two end points of the structure. What are the boundary conditions to be applied on the edges?
Another issue I'm facing is when computing the charge densities by including a non-parabolic coefficient $\alpha$ (for non-parabolic bands), I have to compute these integrals in every iteration - \begin{equation} n = \int^{\infty}_{Ec}F_{1/2}(E)g_n(E)dE \hspace{10pt}p= \int^{E_v}_{-\infty}F_{1/2}(E)g_p(E)dE \end{equation} where $F_{1/2}$ is the 1/2nd order Fermi-distribution and $g(E)$ is the density of states. I have to compute the integrals with $E_c\text{ and }E_v$, but those are unknown. How do I initiate and update their values? How are each of these energies ($E_c,E_v,E_f$) related to each other and the potential of the device?