These are the discretized drift diffusion equations as taken from the book "Analysis and Simulation of Semiconductor Devices" by Siegfried Selberherr. The electron continuity equation is: $$\begin{align} &((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^n.B((\varphi_{i+1,j}-\varphi_{i,j})/V_T).n_{i+1,j} \\ +&((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}^x).D_{i-1/2,j}^n.B((\varphi_{i-1,j}-\varphi_{i,j})/V_T).n_{i-1,j} \\ -&[(((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^n.B((\varphi_{i,j}-\varphi_{i+1,j})/V_T)) \\ +&(((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}).D_{i-1/2,j}^n.B((\varphi_{i,j}-\varphi_{i-1,j})/V_T)) \\ +&(((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^n.B((\varphi_{i,j}-\varphi_{i,j+1})/V_T)) \\ +&(((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^n.B((\varphi_{i,j}-\varphi_{i,j-1})/V_T))].n_{i,j} \\ +&((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^n.B((\varphi_{i,j+1}-\varphi_{i,j})/V_T).n_{i,j+1} \\ +&((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^n.B((\varphi_{i,j-1}-\varphi_{i,j})/V_T).n_{i,j-1} \\ =&R_{i,j}^{EFFECTIVE}.((\Delta_{i-1}^x+\Delta_i^x)/2).((\Delta_{j-1}^y+\Delta_j^y)/2) \end{align}$$
The hole continuity equation is: $$\begin{align} &((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^p.B((\varphi_{i,j}-\varphi_{i+1,j})/V_T).p_{i+1,j} \\ +&((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}^x).D_{i-1/2,j}^p.B((\varphi_{i,j}-\varphi_{i-1,j})/V_T).p_{i-1,j} \\ -&[(((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^p.B((\varphi_{i+1,j}-\varphi_{i,j})/V_T)) \\ +&(((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}^x).D_{i-1/2,j}^p.B((\varphi_{i-1,j}-\varphi_{i,j})/V_T)) \\ +&(((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^p.B((\varphi_{i,j+1}-\varphi_{i,j})/V_T)) \\ +&(((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^p.B((\varphi_{i,j-1}-\varphi_{i,j})/V_T))].p_{i,j} \\ +&((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^p.B((\varphi_{i,j}-\varphi_{i,j+1})/V_T).p_{i,j+1} \\ +&((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^p.B((\varphi_{i,j}-\varphi_{i,j-1})/V_T).p_{i,j-1} \\ =&R_{i,j}^{EFFECTIVE}.((\Delta_{i-1}^x+\Delta_i^x)/2).((\Delta_{j-1}^y+\Delta_j^y)/2 \end{align}$$ where $\Delta_i$ and $\Delta_j$ are the mesh spacing for finite difference grid along X and Y directions. $D_n$,$D_p$ are electron and hole diffusion coefficients.$B(x) = x/(e^x-1)$ is Bernoulli's function and $V_T$ is the thermal voltage.
Also, $$R^{EFFECTIVE} = R^{SRH} + R^{AU} + R^{OPT} + R^{II}$$ where $R^{EFFECTIVE}$ is the effective generation/recombination term. $R^{SRH}$, $R^{AU}$, $R^{OPT}$ and $R^{II}$ are SRH, Auger, Radiative and Impact Ionization Recombination terms.
Now, I am applying the coupled Newton's iterative scheme to solve the above equations where I need to differentiate the equations. When differentiating with respect to potential '$\varphi$', I noticed that $\varphi$ appears in three places namely,
High field Electron/Hole Mobility, $\mu_{n,p}^{high}=\mu_{n,p}^{low}⁄ [1 + (\mu_{n,p}^{LICN}.E_{n,p}/v_{n,p}^{sat})^{\beta_{n,p}}]^{1/\beta_{n,p}}$ and $D_{n,p} = \mu_{n,p}.V_T$;
Bernoulli's function, $B(x)$ where 'x' here means difference between potentials at neighboring grid points and
Impact Ionization $R^{II} = -\alpha_n.(|J_{i,j}^n|/q) - \alpha_p.(|J_{i,j}^p|/q)$, $\alpha_n = a_n.e^{-b_n/|E_{i,j}|}$, $\alpha_p = a_p.e^{-b_p/|E_{i,j}|}$ where $J^{n,p}$ and $E$ represent electron/hole currents and electric field respectively. Current terms are made of potential terms as in, $J_{i+1/2,j}^{nx}=D_{i+1/2,j}^n.[B(\frac{\varphi_{i,j}-\varphi_{i+1,j}}{V_T}).n_{i,j} - B(\frac{\varphi_{i+1,j}-\varphi_{i,j}}{V_T}).n_{i+1,j}]/\Delta_i^x \rightarrow (31)$, $J_{i,j+1/2}^{ny}=D_{i,j+1/2}^n.[B(\frac{\varphi_{i,j}-\varphi_{i,j+1}}{V_T}).n_{i,j} - B(\frac{\varphi_{i,j+1}-\varphi_{i,j}}{V_T}).n_{i,j+1}]/\Delta_j^y \rightarrow (32)$, $J_{i+1/2,j}^{px}=D_{i+1/2,j}^p.[B(\frac{\varphi_{i,j}-\varphi_{i+1,j}}{V_T}).p_{i+1,j} - B(\frac{\varphi_{i+1,j}-\varphi_{i,j}}{V_T}).p_{i,j}]/\Delta_i^x \rightarrow (33)$, $J_{i,j+1/2}^{py}=D_{i,j+1/2}^p.[B(\frac{\varphi_{i,j}-\varphi_{i+1,j}}{V_T}).p_{i,j+1} - B(\frac{\varphi_{i+1,j}-\varphi_{i,j}}{V_T}).p_{i,j}]/\Delta_j^y \rightarrow (34)$
So, basically I need to differentiate all of these to find the Jacobi with respect to $\varphi$. But, the problem arises when the difference between potentials become zero. This means that the entire electron/hole equation becomes a constant when the potential difference is zero. Hence, Jacobi is also zero.
This would mean that the final Jacobian matrix would have entire columns and rows made out of zeros thus making the matrix singular and it is not possible to solve a singular matrix.
Is this the correct approach to solving using the Newton's method. What should I do when the potential difference is zero ?
