Organization

We need organization.

1D Problem

Our one spatial dimension {a≤x≤b} problem is easy to organize:
  Physical:		i=0,1,…,nx   with (Subscript[x, i])=a+i*Δx where Δx=(b-a)/nx
  Linear Algebra:	k=1,…,nx+1 with
  Connection:		k=i+1   
  			or  i=k-1
  Interior Pts:		i=1,…,nx-1
  			or k=2,…,nx
  Boundary Pts:	i=0 (left) and i=nx (right)
  			or k=1(left)  and k=nx+1(right)
  Ghost Pts:		i=-1(left)  and i=nx+1(right)
  (if needed)		or k=0 (left) and k=nx+2(right)
  
  Interior Eq:		1≤i≤nx-1 
  			or 2≤k≤nx
  (explicit)		(Subscript[u, i]^(n+1))=(Subscript[u, i]^n)+νx((Subscript[u, i-1]^n)-2 (Subscript[u, i]^n)+(Subscript[u, i+1]^n))+Δt (Subscript[f, i]^n)
  (implicit)		(Subscript[u, i]^(n+1))=(Subscript[u, i]^n)+νx((Subscript[u, i-1]^(n+1))-2 (Subscript[u, i]^(n+1))+(Subscript[u, i+1]^(n+1)))+Δt (Subscript[f, i]^(n+1))
  (θ)			(Subscript[u, i]^(n+1))=(Subscript[u, i]^n)+ νx (θ((Subscript[u, i-1]^n)-2 (Subscript[u, i]^n)+(Subscript[u, i+1]^n))+(1-θ)((Subscript[u, i-1]^(n+1))-2 (Subscript[u, i]^(n+1))+(Subscript[u, i+1]^(n+1))))+Δt (Subscript[f, i]^(n+θ))
  (Crank-Nicholson)	(Subscript[u, i]^(n+1))=(Subscript[u, i]^n)+ νx ((1/2)((Subscript[u, i-1]^n)-2 (Subscript[u, i]^n)+(Subscript[u, i+1]^n))+(1/2)θ((Subscript[u, i-1]^(n+1))-2 (Subscript[u, i]^(n+1))+(Subscript[u, i+1]^(n+1))))Δt (Subscript[f, i]^(n+1/2))
  
   Boundary Eq:	i=0,nx
   			or k=1,nx+1	
  
  (Dirichlet)		(Subscript[u, 0]^(n+1))=(Subscript[Data, 0]^(n+1))  and/or (Subscript[u, nx]^(n+1))=(Subscript[Data, nx]^(n+1))
  u=Data
  
  (Neumann)		Interior Eq at i=0 with the ghost point eliminated using
  (Subscript[u, x])=NData		 (((Subscript[u, 0]^m)-(Subscript[u, -1]^m))/(Δx))=(Subscript[NData, 0]^m) for m=n,n+1, or n+1/2
  			(((Subscript[u, nx+1]^m)-(Subscript[u, nx]^m))/(Δx))=(Subscript[NData, nx]^m) for m=n,n+1, or n+1/2
  
  (Robin)		Interior Eq at i=0 with the ghost point eliminated using
  ux+αu=RData	 (((Subscript[u, 0]^m)-(Subscript[u, -1]^m))/(Δx))+α (Subscript[u, 0]^m)=(Subscript[RData, 0]^m) for m=n,n+1, or n+1/2
  			 (((Subscript[u, nx+1]^m)-(Subscript[u, nx]^m))/(Δx))+α (Subscript[u, nx]^m)=(Subscript[RData, nx]^m) for m=n,n+1, or n+1/2
  			 
  (Periodic)		Interior Eq at i=0 with the ghost points eliminated using warp around.

To build the matrices we can use the linear algebra index k or the physical index i. The only difference is the start and stop points. When specifying a vector or matrix entry I need to use the linear algebra counting i.  There are several ways to make the matrices for the 1D problem.

Messy (something to plan from)

nx = 4;
AInt=Table[0,{nx+1},{nx+1}];
Do[
k=i+1;
AInt[[k,k-1]]=1.0; AInt[[k,k]]=-2.0; AInt[[k,k+1]]=1.0,
{i,1,nx-1}
]
MatrixForm[AInt]

nx = 21;
{a,b}={0.0,3.14159};
u0[x_]= Sin[2x];
f[x_,t_]= Sin[3x] Sin[t];
κ=0.1;
Δx=((b-a)/nx);
Δt=0.1;
νx= κ((Δt)/((Δx)^2))
AInt=SparseArray[{},{nx+1,nx+1}];
F=Table[0,{nx+1}];
U0=Table[0,{nx+1}];
Do[
k=i+1; (* physical to LA index*)
AInt[[k,k-1]]=1.0; AInt[[k,k]]=-2.0; AInt[[k,k+1]]=1.0;U0[[k]]=u0[a + i Δx], (* Setting IC  *)
{i,1,nx-1}
]
Id=IdentityMatrix[nx+1]; Id[[1,1]]=0; Id[[nx+1,nx+1]]=0;

NIter=33;
Show[
ListLinePlot[U0,
DataRange->{a,b},
PlotRange->{-1,3}],
U=U0;
Do[
F = Table[f[a+i*Δx,(n-1) Δt],{i,0, nx}];
U =(Id + νx AInt).U + Δt*F,
{n,1,NIter}];
ListLinePlot[U,
DataRange->{a,b},
PlotStyle->Red],
PlotLabel->Δt*NIter
]

NIter=35;
Show[
ListLinePlot[U0,
DataRange->{a,b},
PlotRange->{-1,3}],
U=U0;
Do[
F = Table[f[a+i*Δx,(n-1) Δt],{i,0, nx}];
U =(Id + νx AInt).U + Δt*F,
{n,1,NIter}];
ListLinePlot[U,
DataRange->{a,b},
PlotStyle->Red],
PlotLabel->Δt*NIter
]

nx = 21;
{a,b}={0.0,3.14159};
u0[x_]= Sin[2x]+1 + 2 (x/π);
f[x_,t_]= 0 Cos[3x] Sin[t];
κ=0.1;
Δx=((b-a)/nx);
Δt=0.1;
νx= κ((Δt)/((Δx)^2))
AInt=SparseArray[{},{nx+1,nx+1}];
F=Table[0,{nx+1}];
U0=Table[0,{nx+1}];
Do[
k=i+1; (* physical to LA index*)
AInt[[k,k-1]]=1.0; AInt[[k,k]]=-2.0; AInt[[k,k+1]]=1.0;U0[[k]]=u0[a + i Δx], (* Setting IC  *)
{i,1,nx-1}
]
U0[[1]]=u0[a];
U0[[nx+1]]=u0[b];
Id=IdentityMatrix[nx+1]; Id[[1,1]]=0; Id[[nx+1,nx+1]]=0;

NIter=350;
Show[
ListLinePlot[U0,
DataRange->{a,b},
PlotRange->{-1,3}],
U=U0;
Do[
F =Δt* Table[f[a+i*Δx,(n-1) Δt],{i,0, nx}];
F[[1]]=1.0;F[[nx+1]]=3.0;
U =(Id + νx AInt).U + F,
{n,1,NIter}];
ListLinePlot[U,
DataRange->{a,b},
PlotStyle->Red],
PlotLabel->Δt*NIter
]

F[[1]]
F[[1]]

Organized

So I think I will make a couple of functions with obvious names that make rows of the matrix.  I need to make a function for the equations at internal points, the equations at Robin (Neumann is a special case of Robin) boundary points, and periodic points.

For right now I will assume I have the following defined:
	{a, b}, nx, {Δx,Δt}, νx, f, g,α and
	i is the physical spatial counter and n is the physical temporal counter, while
	k is the Linear Algebra counter.
We want to build the sparse array AExp∈(ℝ^((nx+1)×(nx+1))) and vector b∈(ℝ^(nx+1)) to implement the explicit iteration
	(u^new)=(Id + νx A) (u^old)+b
to approximate the solution of 
	(Subscript[u, t])=κ (Subscript[u, xx])+f
of a≤x≤b and t>0 with the intial condition
	u[x,0]=u0[x];
and either periodic conditions 
(P)	u[a-,t]=u[b+,t]
or one of the following boundary conditions at x=a
(D)	u[a, t] = g[a,t] 
(R/N)	(Subscript[u, x])[a, t]+α[a,t] u[a,t] = g[a,t] 
and one of the  boundary conditions at x=ab
(D)	u[b, t] = g[a,t] 
(R/N)	(Subscript[u, x])[b, t]+αb[b,t]u[b,t] = g[b,t].
Note:  α=0 reduces the Robin conditions to a Neumann condition.

We assume that the matrix A and vectors b and U0 are initialied to zero and add the appropriate entries.

InternalEq:= (
k=i+1;
A[[k,k-1]]=  1.0;
A[[k,k      ]]=-2.0;
A[[k,k+1]]=   1.0;
)
InternalRHS:=(
k=i+1;
b[[k]]=f[a+i*Δx,n*Δt]
)
RobinBoundaryEq:=(
k=i+1;
Switch[ i,
0,   mult=  1; x=a,
nx, mult=-1; x=b
,_, Print["Error"]
];
A[[k,k]]=-2.0+(1 + mult α Δx);
A[[k,k+mult]]= 1.0
)
RobinRHS:=(
k=i+1;
Switch[ i,
0,   mult=   1; x=a;
nx, mult=-1; x =b,
_, Print["Error"]
];
b[[k]]=f[x,n*Δt]- mult  Δx g[x,n*Δt];
)
PeriodicBoundaryEq:=(
k=i+1;
Switch[ i,
0, x=a,
nx,x=b,
_, Print["Error"]
];
A[[k,k]]=-2.0;
A[[k,2]]=  1.0;
A[[k,nx-1]]=1.0;
)
PeriodicRHS:=(
k=i+1;
Switch[ i,
0,   mult=   1; x=a,
nx, mult=-1; x =b,
_, Print["Error"]
];
b[[k]]=f[x,n*Δt];
)
IC:=(
k=i+1;
U0[[k]]=u0[a+i*Δx];
)

The stuff up above makes it much cleaner to set up and run a problem

Periodic Boundary Conditions

nx=12;
{a,b}={0.0,3.14159};
{Δx,Δt}={((b-a)/nx),0.01};
κ=0.5;
f[x_,t_]= Sin[x] Cos[t];
νx=κ((Δt)/((Δx)^2))
(* End of Domain intial conditon and forcing function specification *)
(* Intializations *)
n=0;
A= SparseArray[{},{nx+1,nx+1}];
Id = IdentityMatrix[nx+1];
b = Table[0,{nx+1}];
U=b; U0=b;
(* Building matrices and RHS *)
Do[
InternalEq; InternalRHS; IC,
{i,1,nx-1,1}];
Do[
PeriodicBoundaryEq; PeriodicRHS; IC,
{i,0,nx,nx}]
(* Start Run *)
U=U0;
nMax=34;
Do[
U=(Id + νx A).U + b; 
Do[InternalRHS,{i,1,nx-1}]; Do[PeriodicRHS,{i,0,nx,nx}],
{n,1,nMax}]

Dimensions[Id + νx A]
Dimensions[b]

Plotting Matrices and Resolution

{{a,b},{c,d}}={{0,5},{1,4}};
{Δx,Δy}={0.1,0.1};
U = Table[
(x-3)^2 +( y-2)^2,
{x, a, b, Δx},
{y, c, d, Δy}
];
{n,m}=Dimensions[U]
ListPlot3D[U]
ListContourPlot[U]

All our software solves equations A.x=b so we need vectors of unknowns. Easy but confusing :)

u=Flatten[U];
Dimensions[U]
Dimensions[u]
ListLinePlot[u]

UU=Partition[u,31 ];
ListPlot3D[UU]

Counting in a Flattened Array

we need to be able to go backwards and forwards between the Array indexing  and the flat vector indexing.

One way is given {i,j} find k.

Clear[i,j,nx,ny,a,b,c,d,Δx, Δy]
{nx,ny}={2,3};
{{a,b},{c,d}}={{0,1},{2,5}};
Δx = ((b-a)/nx);
Δy = ((d-c)/ny);
x[i_]=a+i*Δx;
y[j_]=c+j*Δy;
k[i_,j_]=i + j*(nx+1)+1;
Show[
Graphics[{
{LightBlue,Rectangle[{a,c},{b,d}]},
Table[{Point[{x[i],y[j]}], 
Text[{i,j}, {x[i],y[j]},{0,1}],
{Red,Text[k[i,j], {x[i],y[j]},{0,-1}]}
},{i,0,nx},{j,0,ny}]
}
],
PlotRange->{{a-Δx,b+Δx},{c-Δy,d+Δy}},
Frame->True,
FrameLabel->{"x","y"},
AspectRatio->1
]

ij[k_]:=Module[
{kk,i,j},
kk = k-1;
i=Mod[kk,nx+1];
j= (kk-i)/(nx+1);
{i,j}
]

Show[
Graphics[{
{LightBlue,Rectangle[{a,c},{b,d}]},
Table[
{i,j}=ij[k];
{Point[{x[i],y[j]}], 
Text[{i,j}, {x[i],y[j]},{0,1}],
{Red,Text[k, {x[i],y[j]},{0,-1}]}
},{k,1,(nx+1)(ny+1)}]
}
],
PlotRange->{{a-Δx,b+Δx},{c-Δy,d+Δy}},
Frame->True,
FrameLabel->{"x","y"},
AspectRatio->1
]

Equations for a 2D Diffusion Equation. 
Periodic and Neumann Conditions

Consider (Subscript[u, t])=k ((Subscript[u, xx])+(Subscript[u, yy]))+f for a≤x≤b, c≤y≤d (with Δx=((b-a)/nx)and Δy=((d-c)/ny)) and initial conditions u(x,y,0)=g(x,y) and boundary conditions u(a,y,t)=u(b,y,t) and (Subscript[u, y])(x,c,t)=(Subscript[u, y])(x,d,t)=0.  Use an explicit scheme with a forward difference in time and centered differences in the spatial variables to approximate the solution at positive times (Subscript[t, n])=n Δt and spatial coordinates (Subscript[x, i])=a+i Δx and (Subscript[y, j])=c+j Δy.

The interior equation at (Subscript[𝒰, i,j]^n)≈u((Subscript[x, i]),(Subscript[y, j]),(Subscript[t, n])) is
	(((Subscript[𝒰, i,j]^(n+1))-(Subscript[𝒰, i,j]^n))/(Δt))=k((((Subscript[𝒰, i-1,j]^n) - 2 (Subscript[𝒰, i,j]^n) + (Subscript[𝒰, i+1,j]^n))/((Δx)^2))+(((Subscript[𝒰, i,j-1]^n) - 2 (Subscript[𝒰, i,j]^n) + (Subscript[𝒰, i,j+1]^n))/((Δy)^2)))+(Subscript[ℱ, i,j]^n)
where (Subscript[ℱ, i,j]^n)≈f((Subscript[x, i]),(Subscript[y, j]),(Subscript[t, n])).

The periodic equation at (Subscript[𝒰, 0,j]^n)≈u((Subscript[x, 0]),(Subscript[y, j]),(Subscript[t, n])) is 
	(((Subscript[𝒰, 0,j]^(n+1))-(Subscript[𝒰, 0,j]^n))/(Δt))=k((((Subscript[𝒰, -1,j]^n) - 2 (Subscript[𝒰, 0,j]^n) + (Subscript[𝒰, 1,j]^n))/((Δx)^2))+(((Subscript[𝒰, 0,j-1]^n) - 2 (Subscript[𝒰, 0,j]^n) + (Subscript[𝒰, 0,j+1]^n))/((Δy)^2)))+(Subscript[ℱ, 0,j]^n)
where we replace (Subscript[𝒰, -1,j]^n)  with (Subscript[𝒰, nx-1,j]^n).

The periodic equation at (Subscript[𝒰, nx,j]^n)≈u((Subscript[x, nx]),(Subscript[y, j]),(Subscript[t, n])) is 
	(((Subscript[𝒰, nx,j]^(n+1))-(Subscript[𝒰, nx,j]^n))/(Δt))=k((((Subscript[𝒰, nx-1,j]^n) - 2 (Subscript[𝒰, nx,j]^n) + (Subscript[𝒰, nx+1,j]^n))/((Δx)^2))+(((Subscript[𝒰, nx,j-1]^n) - 2 (Subscript[𝒰, nx,j]^n) + (Subscript[𝒰, nx,j+1]^n))/((Δy)^2)))+(Subscript[ℱ, nx,j]^n)
where we replace (Subscript[𝒰, nx+1,j]^n)  with (Subscript[𝒰, 1,j]^n).

The equation at (Subscript[𝒰, i,0]^n)≈u((Subscript[x, i]),(Subscript[y, 0]),(Subscript[t, n])) is 
	(((Subscript[𝒰, i,0]^(n+1))-(Subscript[𝒰, i,0]^n))/(Δt))=k((((Subscript[𝒰, i-1,0]^n) - 2 (Subscript[𝒰, i,0]^n) + (Subscript[𝒰, i+1,0]^n))/((Δx)^2))+(((Subscript[𝒰, i,-1]^n) - 2 (Subscript[𝒰, i,0]^n) + (Subscript[𝒰, i,1]^n))/((Δy)^2)))+(Subscript[ℱ, i,0]^n)
where we solve the approximate Neumann condition 
	(((Subscript[𝒰, i,-1]^n) -  (Subscript[𝒰, i,0]^n) )/(Δy))=0
to get 
	(Subscript[𝒰, i,-1]^n) =  (Subscript[𝒰, i,0]^n)
for the ghost value and replace (Subscript[𝒰, i,-1]^n)  with (Subscript[𝒰, 1,0]^n).

The Neumann condition at (Subscript[𝒰, i,ny]^n)≈u((Subscript[x, i]),(Subscript[y, ny]),(Subscript[t, n])) is similar.

There are four corner points.  The equation at (Subscript[𝒰, 0,0]^n)≈u((Subscript[x, 0]),(Subscript[y, 0]),(Subscript[t, n])) should be
	(((Subscript[𝒰, 0,0]^(n+1))-(Subscript[𝒰, 0,0]^n))/(Δt))=k((((Subscript[𝒰, -1,0]^n) - 2 (Subscript[𝒰, 0,0]^n) + (Subscript[𝒰, 1,0]^n))/((Δx)^2))+(((Subscript[𝒰, 0,-1]^n) - 2 (Subscript[𝒰, 0,0]^n) + (Subscript[𝒰, 0,1]^n))/((Δy)^2)))+(Subscript[ℱ, 0,0]^n)
and we can replace the pink values with the Neumann and periodic conditions.

Now the annoying thing is that we somehow want to count from 0 for our discretization and 1 for our matrices etc.

We need to work out how to put these equations into our matrix form
	((OverVector[u])^(n+1))=A·((OverVector[u])^n)+((OverVector[f])^n) 
where (OverVector[u]) and (OverVector[f]) are the flattened versions of the arrays (𝒰^n) and (ℱ^n).

Matrix

The generic internal equation is straightforward.

InternalEq[{i_,j_},{νx_,νy_},{nx_,ny_},k_]:= {
{k[i,j],k[i,j]}->1-2 νx-2νy,
{k[i,j],k[i-1,j]}      ->νx,
{k[i,j],k[i+1,j]}      ->νx,
{k[i,j],k[i,      j-1]}->νy,
{k[i,j],k[i,      j+1]}->νy
}

The periodic equation simply replaces the index of the edge of the domain by the equivalent point on the other edge of the domain. This is easy to code using modular arithmetic.

PeriodicEq[{i_,j_},{νx_,νy_},{nx_,ny_},k_]:={
{k[i,j],k[i,j]}->1-2 νx-2νy,
{k[i,j],k[Mod[i-1,nx],j]}      ->νx,
{k[i,j],k[Mod[i+1,nx],j]}      ->νx,
{k[i,j],k[i,      Mod[j-1,ny]]}->νy,
{k[i,j],k[i,      Mod[j+1,ny]]}->νy
}

The Neumann condition uses the generic homogeneous equation (which introduces a ghost point) and then eliminates the ghost point using the Neumann condition. Right now I do not see how to do all four of the conditions in the same way.

uxNeumannEq[{i_,j_},{νx_,νy_},{nx_,ny_},k_]:={
{k[i,j],k[i,j]}->1-νx-2νy,
Evaluate[
Switch[i,
0,   {k[i,j],k[i+1,j]}->νx,
nx, {k[i,j],k[i-1,j]}->νx]
],
{k[i,j],k[i,      j-1]}->νy,
{k[i,j],k[i,      j+1]}->νy
}
uyNeumannEq[{i_,j_},{νx_,νy_},{nx_,ny_},k_]:=Module[
{},
{
{k[i,j],k[i,j]}->1-2νx-νy,
 {k[i,j],k[i+1,j]}->νx,
{k[i,j],k[i-1,j]}->νx,
Switch[j,
0,  {k[i,j],k[i,j+1]}->νy,
ny,{k[i,j],k[i,j-1]}->νy
]
}
]

The corner equations are periodic in x and Neumann in y  we just need to combine the two things. This looks like this

CornerEq[{i_,j_},{νx_,νy_},{nx_,ny_},k_]:={
{k[i,j],k[i,j]}->1-2νx-νy,
{k[i,j],k[Mod[i+1,nx],j]}->νx,
{k[i,j],k[Mod[i-1,nx],j]}->νx,
Switch[j,
0,  {k[i,j],k[i,j+1]}->νy,
ny,{k[i,j],k[i,j-1]}->νy
]
}

Putting them all together.

A[{νx_,νy_},{nx_,ny_}]:=Module[
{i,j,k,SysSize=(nx+1)*(ny+1)},
k[i_,j_]:= 1+i+j*(nx+1);
SparseArray[
Flatten[
{
Table[  InternalEq[{i,j},{νx,νy},{nx,ny},k], 
 {i,1,nx-1,1},   {j,1,ny-1,1}],
Table[  PeriodicEq[{i,j},{νx,νy},{nx,ny},k], 
 {i,1,nx-1,1},   {j,0,ny,      ny}],
Table[  uyNeumannEq[{i,j},{νx,νy},{nx,ny},k],
  {i,0,nx,       nx},{j,1,ny-1,1}]
(*,
Table[   CornerEq[{i,j},{νx,νy},{nx,ny},k],
  {i,0,nx,       nx},{j,0,ny      ,ny}] *)
}
],
{SysSize, SysSize}]
]

Test

Clear[nx,ny,a,b,c,d,Δx, Δy]
{nx,ny}={6,7};
{{a,b},{c,d}}={{-1.0,1.0},{-1.0,1.0}};
Δx = ((b-a)/nx);
Δy = ((d-c)/ny);
A[{0.1,0.1},{nx,ny}]

Lets make the vector for the RHS of the equation

Clear[f,F]
f[x_,y_]= 1-x^2-y^2
F=Array[0,{SysSize}];
Do[
F[[k[i,j]]]=f[x[i],y[j]],
{i,0,nx},
{j,0,ny}]
FGrid=Partition[F,nx+1];
ListPlot3D[FGrid]

Lets see what the output looks like!

Δt=0.015; κ=2;
{νx,νy}= κ{((Δt)/((Δx)^2)),((Δt)/((Δy)^2))}
U=Table[0, {SysSize}];
Table[
U=A.U+Δt F;
UGrid=Partition[U,{nx+1}];
ListPlot3D[
UGrid,
PlotRange->{0,20},
PlotLabel->i],
{i,1,50}]

Do[
U=A.U+Δt F,
{i,1,500}]

Min[U]

UGrid=Partition[U,{nx+1}];
ListPlot3D[
UGrid,
PlotRange->All,
PlotLabel->(500+5)*Δt]