Sub={(Subscript[u_, i_]):> (Subscript[u, i-1])-2 (Subscript[u, i])+(Subscript[u, i+1])}
Simplify[((Subscript[u, i])/.Sub)/.Sub ]

Clear[MakeA,n]
MakeA[n_]:= SparseArray[{
Band[{1,1}]->6,
Band[{1,2}]->-4,
Band[{1,3}]->1,
Band[{2,1}]->-4,
Band[{3,1}]->1,
Band[{n-1,1}]->1,
Band[{1,n-1}]->1,
{1,n}->-4,
{n,1}->-4
},
{n,n}]

Lets check!

A=MakeA[8];
MatrixForm[A]

CN like scheme.

The von-Neumann stability analysis for Kyle's CN like scheme 
 (((Subscript[U, i]^(m+1))-2(Subscript[U, i]^m)+(Subscript[U, i]^(m-1)))/((Δt)^2))=(1/(2((Δx)^4)))(((Subscript[U, i-2]^(m+1))-4(Subscript[U, i-1]^(m+1))+6(Subscript[U, i]^(m+1))-4(Subscript[U, i+1]^(m+1))+(Subscript[U, i+2]^(m+1)))+((Subscript[U, i-2]^m)-4(Subscript[U, i-1]^m)+6(Subscript[U, i]^m)-4(Subscript[U, i+1]^m)+(Subscript[U, i+2]^m)))
is obtained by substituting (Subscript[U, i]^m)=(λ^m)(E^(2π i (k/n)Δx)) into the up date scheme and solving for λ.

CNScheme=(((Subscript[U, i]^(m+1))-2(Subscript[U, i]^m)+(Subscript[U, i]^(m-1)))/((Δt)^2))==(1/(2((Δx)^4)))(((Subscript[U, i-2]^(m+1))-4(Subscript[U, i-1]^(m+1))+6(Subscript[U, i]^(m+1))-4(Subscript[U, i+1]^(m+1))+(Subscript[U, i+2]^(m+1)))+((Subscript[U, i-2]^m)-4(Subscript[U, i-1]^m)+6(Subscript[U, i]^m)-4(Subscript[U, i+1]^m)+(Subscript[U, i+2]^m)));
StabilityEqn=Simplify[
CNScheme/.(Subscript[U, i_]^(m_)):> (λ^m) (E^(2π I i (k/n)Δx))
]
StabilityEqn=Simplify[λ ((Δt)^2)((StabilityEqn[[1]])/((λ^m) (E^(2π I i (k/n)Δx))))]==Simplify[λ ((Δt)^2) ((StabilityEqn[[2]])/((λ^m) (E^(2π I i (k/n)Δx))))]

StabilityEqn=StabilityEqn/.{
k->((n θ)/(2π  Δx)),
Δt->γ ((Δx)^2)
}

We can rewrite this in the trig form if we want.

ExpToTrig[StabilityEqn]

The question is "How big are the roots of this polynomial as a function of θ and γ?"

G[θ_,γ_]=Abs[λ]/.Solve[StabilityEqn,λ]

Plot3D[
G[θ,γ],
{θ,-π,π},
{γ,0,1},
PlotRange->{0, 1.2}
]

OK, so what do we learn form this!

G[0.1,0.1]

Explicit scheme.

The von-Neumann stability analysis for the explicit scheme 
 (((Subscript[U, i]^(m+1))-2(Subscript[U, i]^m)+(Subscript[U, i]^(m-1)))/((Δt)^2))=(1/((Δx)^4))(((Subscript[U, i-2]^m)-4(Subscript[U, i-1]^m)+6(Subscript[U, i]^m)-4(Subscript[U, i+1]^m)+(Subscript[U, i+2]^m)))
is obtained by substituting (Subscript[U, i]^m)=(λ^m)(E^(2π i (k/n)Δx)) into the up date scheme and solving for λ.

ExpScheme=(((Subscript[U, i]^(m+1))-2(Subscript[U, i]^m)+(Subscript[U, i]^(m-1)))/((Δt)^2))==(1/((Δx)^4))((Subscript[U, i-2]^m)-4(Subscript[U, i-1]^m)+6(Subscript[U, i]^m)-4(Subscript[U, i+1]^m)+(Subscript[U, i+2]^m));
StabilityEqn=Simplify[
ExpScheme/.(Subscript[U, i_]^(m_)):> (λ^m) (E^(2π I i (k/n)Δx))
]
StabilityEqn=Simplify[λ ((Δt)^2)((StabilityEqn[[1]])/((λ^m) (E^(2π I i (k/n)Δx))))]==Simplify[λ ((Δt)^2) ((StabilityEqn[[2]])/((λ^m) (E^(2π I i (k/n)Δx))))]

StabilityEqn=StabilityEqn/.{
k->((n θ)/(2π  Δx)),
Δt->γ ((Δx)^2)
}

We can rewrite this in the trig form if we want.

ExpToTrig[StabilityEqn]

The question is "How big are the roots of this polynomial as a function of θ and γ?"

G[θ_,γ_]=Abs[λ]/.Solve[StabilityEqn,λ]

Plot3D[
G[θ,γ],
{θ,-π,π},
{γ,0,1},
PlotRange->{0, 1.2}
]

OK, so what do we learn form this!

G[0.1,0.1]

Fancy scheme.

The von-Neumann stability analysis for the fancy scheme (with a parameter α)
 (((Subscript[U, i]^(m+1))-2(Subscript[U, i]^m)+(Subscript[U, i]^(m-1)))/((Δt)^2))=(1/((Δx)^4))((α/2)(((Subscript[U, i-2]^(m+1))-4(Subscript[U, i-1]^(m+1))+6(Subscript[U, i]^(m+1))-4(Subscript[U, i+1]^(m+1))+(Subscript[U, i+2]^(m+1)))+((Subscript[U, i-2]^(m-1))-4(Subscript[U, i-1]^(m-1))+6(Subscript[U, i]^(m-1))-4(Subscript[U, i+1]^(m-1))+(Subscript[U, i+2]^(m-1)))) + (1-α)((Subscript[U, i-2]^m)-4(Subscript[U, i-1]^m)+6(Subscript[U, i]^m)-4(Subscript[U, i+1]^m)+(Subscript[U, i+2]^m)))
is obtained by substituting (Subscript[U, i]^m)=(λ^m)(E^(2π i (k/n)Δx)) into the up date scheme and solving for λ.

Clear[α]
FancyScheme=(((Subscript[U, i]^(m+1))-2(Subscript[U, i]^m)+(Subscript[U, i]^(m-1)))/((Δt)^2))==(1/((Δx)^4))((α/2)(((Subscript[U, i-2]^(m+1))-4(Subscript[U, i-1]^(m+1))+6(Subscript[U, i]^(m+1))-4(Subscript[U, i+1]^(m+1))+(Subscript[U, i+2]^(m+1)))+((Subscript[U, i-2]^(m-1))-4(Subscript[U, i-1]^(m-1))+6(Subscript[U, i]^(m-1))-4(Subscript[U, i+1]^(m-1))+(Subscript[U, i+2]^(m-1)))) + (1-α)((Subscript[U, i-2]^m)-4(Subscript[U, i-1]^m)+6(Subscript[U, i]^m)-4(Subscript[U, i+1]^m)+(Subscript[U, i+2]^m)));
StabilityEqn=Simplify[
FancyScheme/.(Subscript[U, i_]^(m_)):> (λ^m) (E^(2π I i (k/n)Δx))
]
StabilityEqn=Simplify[λ ((Δt)^2)((StabilityEqn[[1]])/((λ^m) (E^(2π I i (k/n)Δx))))]==Simplify[λ ((Δt)^2) ((StabilityEqn[[2]])/((λ^m) (E^(2π I i (k/n)Δx))))]

StabilityEqn=StabilityEqn/.{
k->((n θ)/(2π  Δx)),
Δt->γ ((Δx)^2)
}

We can rewrite this in the trig form if we want.

ExpToTrig[StabilityEqn]

The question is "How big are the roots of this polynomial as a function of θ and γ?"

G[θ_,γ_,α_]=Abs[λ]/.Solve[StabilityEqn,λ]

α=0.1;
Plot3D[
G[θ,γ,α],
{θ,-π,π},
{γ,0,1},
PlotRange->{0, 1.2}
]

OK, so what do we learn form this!