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$$
Best Answer
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.