Why does the finite difference script for solving Poisson equation not work for delta function charge? How to fix it?

I know the charge density at the plates of a capacitor, from which the applied voltage and potential profile is to be calculated. Using the finite difference method in MATLAB script I solved the Poisson's equation that gives the correct result for all kind of continuous charge profiles. I checked constant (step function), linearly increasing/decreasing charge profile, in all cases numerical results agree with analytically calculated potential profile. But when I put the delta-function charge to study the original problem (in figure), it gives spooky results. (With $\rho=13.2$$\muC/cm^2, \epsilon_r=30, x=20nm, I would expect the voltage to be 9.92V,while I am getting 40nV!) In the finite different approach: d^2V/dx^2=-\rho(x)/\epsilon in matrix form is Au=h^2$$f$, where $f$=-[$\rho(i)$/$\epsilon$]' and $i$=1,2,......n are the discrete nodes. In my problem, I have $-\rho(1)$=$\rho(n)$=$13.2$$\muC/cm^2, and \rho(i)=0 elsewhere. Therefore I take f=-\rho(i)/\epsilon[-1 0 0 0 .......... 0 0 0 1]'. Then potential is obtained by u=A^{-1}$$f$

My question is, what is the physical reason for such anomaly?

How can I solve this issue and make this code work?  Code:

clc;

clear all;

t=20;  %Thickness value in nm

epsil_r=30;

epsil=epsil_r*8.854*10^-14; %In CGS unit

n=500;              %Number of grid points

x1=0; x2=t*10^-7; %Thickness in cm

x=linspace(x1,x2,n);h=x(2)-x(1);

rho=13.2*10^-6;     %Strength of delta-function charge, C/cm^2

A=diag(-2*ones(1,n))+diag(ones(1,n-1),1)+diag(ones(1,n-1),-1);

A(n,n)=-1;        %Neumann boundary at right node

f=-rho/(epsil)*[-1;zeros(n-2,1);1]; %Right hand side of Poisson.

u=inv(A)*[h^2*f];              %Potential
• Would Computational Science be a better home for this question? – Qmechanic Apr 30 '17 at 16:40
• I think people will be hesitant to discuss the physical behaviour of a simulation whose code they could not verify ... ;) – Sanya Apr 30 '17 at 18:59
• "What is the physical reason for such anomaly?" Either your computer is broken, or (more likely) there is a bug in your code. If you want the code debugged without showing it to anybody you need a mind-reading forum, not SE. – alephzero Apr 30 '17 at 23:48
• Dear altruist, I added the code in the question. I was wondering if there is anything fundamental for which delta-function charge can not be used, or it is all about my code. I am expecting your guideline. – Alam May 2 '17 at 9:32
• I don't understand the code in this question. There's a bunch of definitions, but I don't see what actually acts on those objects. How you plot things is not important, what's important is how you solve things. – Emilio Pisanty May 2 '17 at 10:10

Your units don't add up. You're solving for the equation $$\frac{\mathrm d^2 u}{\mathrm dx^2} = f$$ where $f = \rho/\epsilon$. If you take your solution at face value, then $u$ will have the dimensionality \begin{align} [u] & = [h^2 f] = [h^2\rho/\epsilon_0] = [Q/\epsilon_0]. \end{align} Since $[\epsilon_0]=[Q/(VL)]$, this means that $[u]=[VL]$, a voltage times length, which is obviously not what you intended. Since the units don't match up, it is meaningless to talk about inconsistencies in the numbers - if you did this in some other unit system you would get a different answer.
The problem is in your handling of the surface charge associated with the delta function. The right-hand side of the Poisson equation requires a legitimate volumetric charge density, or something dimensionally equivalent to that. If $\sigma$ is a surface charge density, with $[\sigma]=[Q/L^2]$, then $$\rho(x,y,z) = \sigma\delta(x-x_0)$$ is a correct volumetric charge density because $[\delta(x-x_0)] = [1/L]$, but your discretization procedure has dropped this fact out.
To fix this, I would suggest dropping the pretense that you have a true singular charge density and instead model it as a finite slab of thickness $d$ and volumetric charge density $\sigma/d$, where you can set $d$ to your discretization length if you want. Ultimately, the value of $d$ should be unimportant to your finite result - but it's crucial that you get your units right.