Cranck Nicolson 1 and 2 Dimensions

The plan is to implement the CN scheme is 1D for

Differencing Matrices in 1D

The differencing (Subscript[δ, 2]) matrix for differencing with peridic boundary conditions and nx internal points and one boundary point.

DiffMat1DPeriodic[nx_] := SparseArray[{
Band[{1,1}]->-2.0,
Band[{1,2}]->1.0,
Band[{2,1}]->1.0,
{1,nx+1}->1.0,
{nx+1,1}->1.0
},
{nx+1, nx+1}]
MatrixForm[DiffMat1DPeriodic[6]]

Now the differencing matrix for Dirichlet conditions with nx internal points and two boundary point.

DiffMat1DDirichlet[nx_] := SparseArray[{
Band[{2,2},{nx+1,nx+2}]->-2.0,
Band[{2,3},{nx+1,nx+2}]->1.0,
Band[{2,1},{nx+1,nx+2}]->1.0
},
{nx+2, nx+2}]
MatrixForm[DiffMat1DDirichlet[6]]

Time Steps in 1D

The plan is to advance nt time steps

CNStepPeriodic[UIn_,nt_,ν_]:= Module[
{nx=Length[UIn]-1,Ap, Am,U=UIn},
Ap= IdentityMatrix[nx+1]+0.5*ν DiffMat1DPeriodic[nx+1];
Am= IdentityMatrix[nx+1]-0.5*ν DiffMat1DPeriodic[nx+1];
Do[
U=LinearSolve[Ap,Am.U],
{nt}];
U
]

The differencing (Subscript[δ, 2]) matrix for differencing with peridic boundary conditions and n internal points and one boundary point.

Code

Clear[Relax3,ω,νx,νy,νz,U,F]
SetAttributes[Relax3,HoldRest]
Relax3 [{ω_,{νx_,νy_,νz_}},U_,F_]:=Module[
{nx,ny,nz,i,j,k,α,β},
{nx,ny,nz}=Dimensions[F];
α=0.5ω/(νx+νy+νz);
β = 1-ω;
(* red cells i+j+k even first *) 
Do[
U[[i,j,k]]=α (F[[i,j,k]] 
+ νx (U[[i-1,j,k]]+U[[i+1,j,k]])
+ νy (U[[i,j-1,k]]+U[[i,j+1,k]])
+ νz (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + β U[[i,j,k]],
{i,2,nx-1,2}, (* even *)
{j,2,ny-1,2}, (* even *)
{k,2,nz-1,2}     (* even *)
];
Do[
U[[i,j,k]]=α (F[[i,j,k]] 
+ νx (U[[i-1,j,k]]+U[[i+1,j,k]])
+ νy (U[[i,j-1,k]]+U[[i,j+1,k]])
+ νz (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + β U[[i,j,k]],
{i,2,nx-1,2},  (* even *)
{j,3,ny-1,2},  (* odd *)
{k,3,nz-1,2}     (* odd *)
];
Do[
U[[i,j,k]]=α (F[[i,j,k]] 
+ νx (U[[i-1,j,k]]+U[[i+1,j,k]])
+ νy (U[[i,j-1,k]]+U[[i,j+1,k]])
+ νz (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + β U[[i,j,k]],
{i,3,nx-1,2},  (* odd *)
{j,2,ny-1,2},  (* even *)
{k,3,nz-1,2}     (* odd *)
];
Do[
U[[i,j,k]]=α (F[[i,j,k]] 
+ νx (U[[i-1,j,k]]+U[[i+1,j,k]])
+ νy (U[[i,j-1,k]]+U[[i,j+1,k]])
+ νz (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + β U[[i,j,k]],
{i,3,nx-1,2},  (* odd *)
{j,3,ny-1,2},  (* odd *)
{k,2,nz-1,2}     (* even *)
];
(* Now black Cells i+j+k odd *)
Do[
U[[i,j,k]]=α (F[[i,j,k]] 
+ νx (U[[i-1,j,k]]+U[[i+1,j,k]])
+ νy (U[[i,j-1,k]]+U[[i,j+1,k]])
+ νz (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + β U[[i,j,k]],
{i,3,nx-1,2}, (* odd *)
{j,2,ny-1,2}, (* even *)
{k,2,nz-1,2}     (* even *)
];
Do[
U[[i,j,k]]=α (F[[i,j,k]] 
+ νx (U[[i-1,j,k]]+U[[i+1,j,k]])
+ νy (U[[i,j-1,k]]+U[[i,j+1,k]])
+ νz (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + β U[[i,j,k]],
{i,2,nx-1,2}, (* even *)
{j,3,ny-1,2}, (* odd *)
{k,2,nz-1,2}     (* even *)
];
Do[
U[[i,j,k]]=α (F[[i,j,k]] 
+ νx (U[[i-1,j,k]]+U[[i+1,j,k]])
+ νy (U[[i,j-1,k]]+U[[i,j+1,k]])
+ νz (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + β U[[i,j,k]],
{i,2,nx-1,2}, (* even *)
{j,2,ny-1,2}, (* even *)
{k,3,nz-1,2}     (* odd *)
];
Do[
U[[i,j,k]]=α (F[[i,j,k]] 
+ νx (U[[i-1,j,k]]+U[[i+1,j,k]])
+ νy (U[[i,j-1,k]]+U[[i,j+1,k]])
+ νz (U[[i,j,k-1]]+U[[i,j,k+1]])
 ) + β U[[i,j,k]],
{i,3,nx-1,2}, (* odd *)
{j,3,ny-1,2}, (* odd *)
{k,3,nz-1,2}     (* odd *)
];
]

Runs

The 2D equation (with Δx≠Δy) is
	(((Subscript[U, i-1,j])-2(Subscript[U, i,j])+(Subscript[U, i+1,j]))/((Δx)^2))+(((Subscript[U, i,j-1])-2(Subscript[U, i,j])+(Subscript[U, i,j+1]))/((Δy)^2))=(Subscript[F, i,j])
The 3D equation (with Δx≠Δy) is
	(((Subscript[U, i-1,j,k])-2(Subscript[U, i,j,k])+(Subscript[U, i+1,j,k]))/((Δx)^2))+(((Subscript[U, i,j-1,k])-2(Subscript[U, i,j,k])+(Subscript[U, i,j+1,k]))/((Δy)^2))+(((Subscript[U, i,j,k-1])-2(Subscript[U, i,j,k])+(Subscript[U, i,j,k+1]))/((Δz)^2))=(Subscript[F, i,j,k]).

The relaxation iteration equations are derived as for 2D
	νx((Subscript[U, i-1,j])-2(Subscript[U, i,j])+(Subscript[U, i+1,j]))+νy((Subscript[U, i,j-1])-2(Subscript[U, i,j])+(Subscript[U, i,j+1]))=(Subscript[F, i,j])
where νx=((Δx)^(-2)) and νy=((Δy)^(-2)) etc.
	2(νx+νy)(Subscript[U, i,j])=(Subscript[F, i,j])+νx((Subscript[U, i-1,j])+(Subscript[U, i+1,j]))+νy((Subscript[U, i,j-1])+(Subscript[U, i,j+1]))
the non-relaxed (Gauss-Seidel)  iteration is 
	(Subscript[U, i,j])=(1/(2(νx+νy)))((Subscript[F, i,j])+νx((Subscript[U, i-1,j])+(Subscript[U, i+1,j]))+νy((Subscript[U, i,j-1])+(Subscript[U, i,j+1])))
the relaxed iteration is
	(Subscript[U, i,j])=ω((1/(2(νx+νy)))((Subscript[F, i,j])+νx((Subscript[U, i-1,j])+(Subscript[U, i+1,j]))+νy((Subscript[U, i,j-1])+(Subscript[U, i,j+1]))))+(1-ω)(Subscript[U, i,j])
The 3D version is similar
	(Subscript[U, i,j,k])=ω((1/(2(νx+νy+νz)))((Subscript[F, i,j,k])+νx((Subscript[U, i-1,j,k])+(Subscript[U, i+1,j,k]))+νy((Subscript[U, i,j-1,k])+(Subscript[U, i,j+1,k]))+νz((Subscript[U, i,j,k-1])+(Subscript[U, i,j,k+1]))))+(1-ω)(Subscript[U, i,j,k])

All of these are easy to code in place and perform in a red/black ordering

{bx,by,bz}=N[{π,π,π}];
{nx,ny,nz}={32,32,23};
{Δx,Δy,Δz}={bx,by,bz}/{nx+1,ny+1,nz+1}
U=ConstantArray[0,{nx+2,ny+2,nz+2}];
f[x_,y_,z_]:= Sin[x] Sin[y] Sin[z];
F=Table[
f[x,y,z],
{x,0,bx,Δx},
{y,0,by,Δy},
{z,0,bz,Δz}
];
ω=1.8;
{νx,νy,νz}=1/(({Δx,Δy,Δz})^2);
Do[
Relax3[{ω,{νx,νy,νz}},U,F],
{120}]

Err=U-F/3;
Max[Abs[Err]]

ListContourPlot[ U[[6]]]

ListContourPlot3D[U,
Contours->{0.1}]
ListContourPlot3D[F,
Contours->{0.1}]

νx

Relax3[{ω,{νx,νy,νz}},U_,F_]

?Relax3