I want to solve the diffusion equation using the method of lines with Neumann boundary conditions $$ \frac{\partial p}{\partial t}=\frac{\partial^2p}{\partial x^2}\\ \frac{\partial p}{\partial x}(x=0)=\frac{\partial p}{\partial x}(x=1)=0\\ p(x,0)=\begin{cases} 1/dx \quad x=0.5,\\ 0 \quad elsewhere \end{cases} $$ For sake of simplicity I use the forward Euler so I choose $dx=0.1, \; dt=0.001$ to ensure the stability. Using a second order scheme I solve the boundary conditions imposing $$ \frac{p_1-p_{-1}}{2 dx}=0,\\ \frac{p_{N+1}-p_{N-1}}{2 dx}=0.\\ $$ So the tridiagonal matrix looks like $$ \frac{1}{dx^2} \begin{pmatrix} -2&2&0&0\\ 1&-2&1&0\\ 0&1&-2&1\\ 0&0&2&-2 \end{pmatrix}. $$ At this point the solution of my equation is given by $$ \vec{p}^{k+1}=dt(T\vec{p}^k)+\vec{p}^k. $$ This boundary conditions impose the conservation of probability but it seems that is not the case since for different $N$ (number of grid points on the $x$ axis) I obtain different probabilities. For example, for $N=10,20,100,1000$ I obtained that the total probability is $1.11111,1.05263,1.00012,1.0000$. It seems that I need a really fine grid to obtain the correct physical result.
Did I make a conceptual mistake or I just made an error in my code? Thank for the answers