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$$
 A: 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. 
