The app below solves the Poisson equation to determine the band bending, the charge distribution, and the electric field in a MOS capacitor with a ptype substrate.
Band diagram

 
Charge density

Electric field

In this app, the bulk of the semiconductor is grounded and $V$ is the voltage on the metal gate, $t_{\text{ox}}$ is the thickness of the oxide, $E_g$ is the band gap which can be temperature dependent, $\epsilon_{\text{ox}}$ is the relative dielectric constant of the oxide, $\epsilon_s$ is the relative dielectric constant of the semiconductor, $E_s$ is the electric field at the semiconductoroxide interface, $Q$ is the charge on the charge on the semiconductor, $\phi_m$ is the work function of the metal, $\chi_s$ is the electron affinity of the semiconductor, $V_s$ is the voltage at the semiconductoroxide interface, $T$ is the absolute temperature, $N_A$ is the concentration of acceptors, $N_c(300)$ is the effective density of states in the conduction band at 300 K, and $N_v(300)$ is the effective density of states in the valence band at 300 K. Note that in the band diagram, the band gap of the oxide is drawn 2 eV larger than the band gap of the semiconductor. In reality, the band gap of the oxide would be much larger but it is easier to see the band bending when the diagram is plotted with a small band gap for the gate oxide.
The and buttons change the gate voltage by $\pm 0.01$ V.
When the metal and the semicondutor are electrically connected, electrons flow from the material with a low work function to the material with a high work function. This current flow stops when a builtin voltage is established that is equal to the work function difference, $V_{bi}=\phi_m\phi_s$. The work function of a semiconductor depends on the electron affinity $\chi_s$ and the doping, $\phi_s=\chi_s+E_cE_F$.
A MOS Capacitor can be in three regimes accumulation, depletion, and inversion. [1,2] The boundary between accumulation and depletion is the flatband voltage and the boundary between depletion and inversion is the threshold voltage. The flatband voltage is the voltage where no band bending occurs, $V_{fb}=V_{bi}=\phi_m\phi_s$. At the flatband voltage, the bands are flat. Since the electric field is the derivative of the band, the electric field is zero everywhere. The hole concentration $p$ is the same as the acceptor concentration and the charge density $\rho$ is zero. Press the '+' and '' buttons to move into accumulation and depletion.
Accumulation occurs when the gate voltage $V$ is negative so that holes are attracted to the oxide interface. The valence band bends up towards the Fermi energy so that the hole concentration, $p$ increases near the oxide. The black line on the left that represents the Fermi energy in the metal moves up for negative voltages because the band diagram plots the electron energy and the energy of electrons increases as a negative voltage is applied.
In the depletion regime, a positive gate voltage pushes the mobile holes away from the oxide and the negatively charged acceptors are left behind. The valence band bends away from the Fermi energy at the oxide corresponding to a lower hole concentration there. The negative charge in the semiconductor is compensated by a positive charge on the surface of the metal. This is indicated by the arrow in the charge plot. As the gate voltage is made more positive, the depletion width grows and the bands bend down. Eventually, the conduction band gets closer to the Fermi energy that the valence band. The condition that $n > p$ at the oxide is called weak inversion. Strong inversion is the condition that $n = N_A$ at the oxide. Strong inversion occurs at the threshold voltage $V_T$. At this point a layer of mobile electrons appears at the semiconductor/oxide nterface. This is called the inversion channel. For voltages $V > V_T$ the depletion width does not increase any further and any extra positive charge on the surface of the metal is compensated by electrons in the inversion channel. Notice that in inversion, the electric field in the semiconductor stays constant; it only increases in the oxide.
To calculate the band bending, we start with Gauss's law,
\[ \begin{equation} \nabla \cdot \vec{E} = \frac{\rho}{\epsilon_s\epsilon_0}. \end{equation} \]Combining this with $\vec{E}=\nabla V$ yields the Poisson equation,
\[ \begin{equation} \nabla^2V = \frac{\rho}{\epsilon_s\epsilon_0}, \end{equation} \]where, for a MOS capacitor with a ptype substrate, the charge density is $\rho = e\left(N_An+p\right)$ and the charge carrier concentrations are,
\begin{equation} n=N_c(300)\left(\frac{T}{300}\right)^{3/2}\exp\left(\frac{E_FE_c}{k_BT}\right)\qquad \text{and}\qquad p=N_v(300)\left(\frac{T}{300}\right)^{3/2}\exp\left(\frac{E_vE_F}{k_BT}\right). \end{equation}Using the relation $e\frac{dV}{dx} = \frac{E_v}{dx}$ the Poisson equation can be written as a second order differential equation for $E_v(x)$,
\[ \begin{equation} \frac{d^2E_v}{dx^2} = \frac{e^2}{\epsilon_s\epsilon_0}\left(N_AN_c\exp\left(\frac{E_gE_v}{k_BT}\right)+N_v\exp\left(\frac{E_v}{k_BT}\right)\right). \end{equation} \]This differential equation was solved numerically using the shooting method [3]. First the maximum depletion width max$(x_p)$ and the threshold voltage $V_T$ are estimated using the analytic formulas from the depletion approximation.
\[ \begin{equation} x_p = 2\sqrt{\frac{\epsilon_{\text{semi}}\epsilon_0 k_BT}{e^2N_A}\ln\left(\frac{N_A}{n_i}\right)}. \end{equation} \] \begin{equation} V_T = \frac{2t_{ox}}{\epsilon_{ox}}\sqrt{\epsilon_{\text{semi}}N_Ak_BT \ln \left (\frac{N_A}{n_i} \right )} +\frac{2k_BT}{e} \ln \left (\frac{N_A}{n_i} \right ) +V_{fb} \end{equation}Far from the oxide, the valence band satisfies the conditions $E_{v} = k_BT\ln\left(\frac{N_A}{N_v}\right)=E_{v0}$ and $\frac{dE_v}{dx}=0$. To determine the band bending, we start a distance of $1.8x_p$ from the oxide with $E_{v} = k_BT\ln\left(\frac{N_A}{N_v}\right)=E_{v0}$ and a small value of $\frac{dE_v}{dx}=0$. The Poisson equation is integrated numerically using the midpoint method until the semiconductor oxide interface. This gives us the voltage $V_s$ at the semiconductor/oxide interface and the electric field $E_s$ at that point. The voltage on the gate is,
\begin{equation} V = \frac{\epsilon_{\text{semi}}E_s}{\epsilon_{\text{ox}}}t_{\text{ox}}+V_s. \end{equation}This is the correct gate voltage for the boundary conditions we chose on the right but generally will not be the gate voltage we wanted. The position where the integration starts from is then adjusted right or left and the integration is performed again until the integration yields the correct gate voltage. The voltage calculated by numerical integration is printed out as $V_{\text{shoot}}$.
The app gives erroneous results if the valence band or the conduction band get closer than about $3k_BT$ from the Fermi energy because the formulas for $n$ and $p$ are only valid when the valence and conductance bands are not too close to the Fermi energy.