# Neumann boundary conditions in finite difference

Inspired by this question, the finite difference solution for the PDE of $$u_t = \kappa u_{xx}$$ with initial/boundary conditions of $$u(x,0) = 0\\ u(0,t)=100\\ u_x(l,t)=A$$

is $$T[n+1,j] = T[n,j] + s(T[n,j-1] - 2T[n,j] + T[n, j-1] + 2Adx)$$

what is the solution if both boundary conditions are Neumann $$u_x(0,t)=A\\ u_x(l,t)=0$$

Consider the grid $x_{0}<\cdots<x_{M}$ where $x_{j}=j\Delta x$ and $x_{M}=M\Delta x=l$. Discretizing the PDE using implicit-in-time and central-in-space differences yields (for $0<j<M$) $$\frac{u_{j}^{n}-u_{j}^{n-1}}{\Delta t}=\kappa\frac{u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}}{(\Delta x)^{2}}$$ (I am using $u_{j}^{n}\equiv T[n,j]$ to make the notation simpler). Simplifying, $$u_{j}^{n}+\frac{\kappa\Delta t}{(\Delta x)^{2}}\left(-u_{j-1}^{n}+2u_{j}^{n}-u_{j+1}^{n}\right)=u_{j}^{n-1}.$$ As for the boundary conditions ($j=0,M$), we have $$\frac{u_{1}^{n}-u_{0}^{n}}{\Delta x}=A\qquad\text{and}\qquad\frac{u_{M}^{n}-u_{M-1}^{n}}{\Delta x}=B$$ (I have used $B$ instead of directly substituting $B=0$ to keep things general). Equivalently, $$u_{0}^{n}-u_{1}^{n}=-A\Delta x\qquad\text{and}\qquad u_{M}^{n}-u_{M-1}^{n}=B\Delta x.$$ Let $u^{n}=(u_{0}^{n},\ldots,u_{M}^{n})^{\intercal}$ be the vector consisting of the solution at time $n\Delta t$. Then, we can summarize all of the above into the linear system $$\mathcal{M}u^{n}=r^{n}$$ where $\mathcal{M}$ is the matrix $$I+\frac{\kappa\Delta t}{(\Delta x)^{2}}\begin{pmatrix}0 & 0\\ -1 & 2 & -1\\ & -1 & 2 & -1\\ & & \ddots & \ddots & \ddots\\ & & & -1 & 2 & -1\\ & & & & -1 & 2 & -1\\ & & & & & & 0 & 0 \end{pmatrix}-\begin{pmatrix}0 & -1\\ 0 & 0 & 0\\ & 0 & 0 & 0\\ & & \ddots & \ddots & \ddots\\ & & & 0 & 0 & 0\\ & & & & 0 & 0 & 0\\ & & & & & -1 & 0 \end{pmatrix}$$ and $r^{n}=(r_{0}^{n},\ldots,r_{M}^{n})^{\intercal}$ is the vector with entries $$r_{j}^{n}=\begin{cases} -A\Delta x & \text{if }j=0\\ u_{j}^{n-1} & \text{if }0<j<M\\ B\Delta x & \text{if }j=M. \end{cases}$$ It is possible to show that $\mathcal{M}$ is a nonsingular M-matrix (you can use this to do so). Note that you still need to specify a boundary condition at $t=0$, i.e. $$u(\cdot,0)=\varphi(\cdot).$$ Once this is specified, we can set $$u_{j}^{0}=\varphi(x_{j})$$ and compute $u^{1}=\mathcal{M}^{-1}b^{1}$, $u^{2}=\mathcal{M}^{-1}b^{2}$, etc.