(SubPlus[ψ])[x_]:=(Subscript[A, n]) Exp[I k x] + (Subscript[B, n]) Exp[-I k x]
(SubMinus[ψ])[x_]:=(Subscript[A, n+1]) Exp[I k x] +(Subscript[B, n+1]) Exp[-I k x]
(Subscript[V, n_])[y_]:=α*(Subscript[δ, x])[y-(Subscript[x, n])]
(Subscript[X, n_]):=n a
" "
H[t_]:=Solve[{D[(SubMinus[ψ])[(Subscript[x, n])],(Subscript[x, n])]  -  D[(SubPlus[ψ])[(Subscript[x, n])],(Subscript[x, n])]==((2m)/(ℏ^2))α (SubPlus[ψ])[(Subscript[x, n])]&&(SubMinus[ψ])[(Subscript[x, n])]==(SubPlus[ψ])[(Subscript[x, n])], β==((m α)/((ℏ^2)k))},{(Subscript[A, n+1]),(Subscript[B, n+1])},{α,m,ℏ}][[1]]/.n->t

Firstly, let’s see the potential that we are dealing with:

Manipulate[Plot[(∑_(l=1)^n)((1.*^3)/(π*((((s-(Subscript[X, 10l])))^2)+((1.*^-3)^2))))/.a->1,{s,-020,120},PlotRange->Automatic,PlotStyle->{AbsoluteThickness[3.],Black}],{n,0,10,1}]

Starting from T.I. Schrodinger’s Equation

Dt[Ψ[x],{x,2}]==((-2m)/(ℏ^2))(T-(Subscript[V, n])[x])Ψ[x]//TraditionalForm
Integrate[Dt[Ψ[x],{x,2}]==-((2m)/(ℏ^2))(T-(Subscript[V, n])[x])Ψ[x],{x,-ϵ,ϵ}]//TraditionalForm
Dt[Ψ[ϵ],{ϵ,1}]+Dt[Ψ[-ϵ],{ϵ,1}]==((2m)/(ℏ^2))α Ψ[(Subscript[x, n])]//TraditionalForm
Limit[Dt[Ψ[ϵ],{ϵ,1}]+Dt[Ψ[-ϵ],{ϵ,1}],ϵ->(Subscript[x, n])]==(SubPlus[Ψ])[(Subscript[x, n])]-(SubMinus[Ψ])[(Subscript[x, n])]//TraditionalForm
(SubPlus[Ψ])'[(Subscript[x, n])]-(SubMinus[Ψ])'[(Subscript[x, n])]==((2m)/(ℏ^2))α Ψ[(Subscript[x, n])]//TraditionalForm

So we have two boundary conditions, The first is:

(SubPlus[Ψ])[(Subscript[x, n])]==(SubMinus[Ψ])[(Subscript[x, n])]//TraditionalForm
(SubPlus[ψ])[x]==(SubMinus[ψ])[x]//TraditionalForm

The second:

(SubPlus[Ψ])'[(Subscript[x, n])]-(SubMinus[Ψ])'[(Subscript[x, n])]==((2m)/(ℏ^2))Ψ[(Subscript[x, n])]//TraditionalForm
D[(SubMinus[Ψ])[x],x]  -  D[(SubPlus[Ψ])[x],x]==((2m)/(ℏ^2))α ((SubPlus[ψ])[(Subscript[x, n])])/.x->(Subscript[x, n])//TraditionalForm
D[(SubPlus[ψ])[x],x]==D[(SubMinus[ψ])[x],x]==((2m)/(ℏ^2))α ((SubPlus[ψ])[(Subscript[x, n])])/.x->(Subscript[x, n])//TraditionalForm

Calling β = ((m α)/((ℏ^2)k)), Then simplifying the last equation

D[(SubMinus[ψ])[x],x]  -  D[(SubPlus[ψ])[x],x]==((2m)/(ℏ^2))α ((SubPlus[ψ])[(Subscript[x, n])])/.x->(Subscript[x, n])//TraditionalForm//Simplify
  (E^(-I k (Subscript[x, n]))) ((Subscript[A, n]) (-(E^(2 I k (Subscript[x, n]))))+(Subscript[B, n])+F (E^(2 I k (Subscript[x, n])))-G)==((2 α m ((SubPlus[ψ])[(Subscript[x, n])]))/(I k (ℏ^2)))//TraditionalForm
  (E^(-I k (Subscript[x, n]))) ((Subscript[A, n]) (-(E^(2 I k (Subscript[x, n]))))+(Subscript[B, n])+F (E^(2 I k (Subscript[x, n])))-G)==-2I β ((SubPlus[ψ])[(Subscript[x, n])])//TraditionalForm

Now we have these two equations that we will solve for (Subscript[A, n+1]) and (Subscript[B, n+1]) then construct the matrix (Subscript[M, n]) which is the coefficient matrix of (Subscript[A, n]) and (Subscript[B, n]) and represent the Transfer Matrix

(SubMinus[ψ])[(Subscript[x, n])]==(SubPlus[ψ])[(Subscript[x, n])]//TraditionalForm //TraditionalForm
 (E^(-I k (Subscript[x, n]))) ((Subscript[A, n]) (-(E^(2 I k (Subscript[x, n]))))+(Subscript[B, n])+F (E^(2 I k (Subscript[x, n])))-G)==-2I β ((Subscript[A, n])+(Subscript[B, n]))//TraditionalForm

S = Solve[{D[(SubMinus[ψ])[(Subscript[x, n])],(Subscript[x, n])]  -  D[(SubPlus[ψ])[(Subscript[x, n])],(Subscript[x, n])]==((2m)/(ℏ^2))α (SubPlus[ψ])[(Subscript[x, n])]&&(SubMinus[ψ])[(Subscript[x, n])]==(SubPlus[ψ])[(Subscript[x, n])], β==((m α)/((ℏ^2)k))},{(Subscript[A, n+1]),(Subscript[B, n+1])},{α,m,ℏ}][[1]];
S/.Rule->Equal//TraditionalForm
"(Subscript[M, n])"==MatrixForm@Normal@CoefficientArrays[{(Subscript[A, n+1])/. S,(Subscript[B, n+1])/. S},{(Subscript[A, n]),(Subscript[B, n])}][[2]]//TraditionalForm

Arranging the matrices yields

MatrixForm@{(Subscript[A, n+1]),(Subscript[B, n+1])}=="(Subscript[M, n])"MatrixForm@{(Subscript[A, n]),(Subscript[B, n])}//TraditionalForm
MatrixForm@{(Subscript[A, n+1]),(Subscript[B, n+1])} ==MatrixForm@Normal@CoefficientArrays[{(Subscript[A, n+1])/. S,(Subscript[B, n+1])/. S},{(Subscript[A, n]),(Subscript[B, n])}][[2]]| MatrixForm@{(Subscript[A, n]),(Subscript[B, n])}//MatrixForm//TraditionalForm

Let’s Assume That we have three barriers, so n = 0, 1 , 2. We will now observe the behavior of the relationships between coefficients

MatrixForm@{(Subscript[A, n+1]),(Subscript[B, n+1])}=="(Subscript[M, n])"MatrixForm@{(Subscript[A, n]),(Subscript[B, n])}//TraditionalForm
MatrixForm@{(Subscript[A, n+1]),(Subscript[B, n+1])} ==Normal@CoefficientArrays[{(Subscript[A, n+1])/. S,(Subscript[B, n+1])/. S},{(Subscript[A, n]),(Subscript[B, n])}][[2]]| MatrixForm@{(Subscript[A, n]),(Subscript[B, n])}//MatrixForm//TraditionalForm

Setting n = 2, 1 ,0. Respectively, we get:

MatrixForm@{(Subscript[A, n+1]),(Subscript[B, n+1])}==(Subscript[M, n])MatrixForm@{(Subscript[A, n]),(Subscript[B, n])}/.n->2//TraditionalForm
MatrixForm@{(Subscript[A, n+1]),(Subscript[B, n+1])}==(Subscript[M, n])MatrixForm@{(Subscript[A, n]),(Subscript[B, n])}/.n->1//TraditionalForm
MatrixForm@{(Subscript[A, n+1]),(Subscript[B, n+1])}==(Subscript[M, n])MatrixForm@{(Subscript[A, n]),(Subscript[B, n])}/.n->0//TraditionalForm

Combining them in one equation yields:

MatrixForm@{(Subscript[A, n+1]),(Subscript[B, n+1])}==((∏_(l=0)^2)(Subscript[M, l]))MatrixForm@{(Subscript[A, 0]),(Subscript[B, 0])}/.n->2//TraditionalForm

We can now generalize this for n = 0, 1 , 2 ..., N:

MatrixForm@{(Subscript[A, N+1]),(Subscript[B, N+1])}==(∏_(n=0)^N) (Subscript[M, n])MatrixForm@{(Subscript[A, 0]),(Subscript[B, 0])}//TraditionalForm

Here you can Play with N to see the resultant matrix:

Manipulate[MatrixForm@{(Subscript[A, N+1]),(Subscript[B, N+1])}==((∏_(n=0)^N) (Subscript[M, n]))MatrixForm@{(Subscript[A, 0]),(Subscript[B, 0])}//TraditionalForm,{N,0,10,1}]

Now let’s check the determinant of (Subscript[M, n]) and check if it yields unity

"(Subscript[M, n])"==MatrixForm@Normal@CoefficientArrays[{(Subscript[A, n+1])/. H[n],(Subscript[B, n+1])/. H[n]},{(Subscript[A, n]),(Subscript[B, n])}][[2]]//TraditionalForm
M[n_]:=Normal@CoefficientArrays[{(Subscript[A, n+1])/. H[n],(Subscript[B, n+1])/. H[n]},{(Subscript[A, n]),(Subscript[B, n])}][[2]];
De[m_]:=Reverse@Table[M[n],{n,0,m,1}]
wq[qq_]:=ParallelCombine[Dot[Sequence@@##]&,De[qq]/.k->m/.β->(1/m)/.(Subscript[x, n_])->n,Dot]
"M =\!\(\*UnderoverscriptBox[\(\(\\ \)\(∏\)\), \(n = 0\), \(N\)]\) (Subscript[M, n])"//TraditionalForm
( Det[ M]//HoldForm)==Det[wq[1]]//FullSimplify//TraditionalForm

"The determinant is:" Manipulate[Det[wq[N]]//FullSimplify,{N,0,15,1}]

From previous calculations of the B.C. we obtain these relations between the Coefficients:

(Subscript[A, n+1])==-I (β (Subscript[A, n])+I (Subscript[A, n])+β (Subscript[B, n]) (E^(-2 I k (Subscript[x, n]))))//TraditionalForm
(Subscript[B, n+1])==I β (Subscript[A, n]) (E^(2 I k (Subscript[x, n])))+I β (Subscript[B, n])+(Subscript[B, n])//TraditionalForm

Taking two wave functions, one represents far right incedent wave with (Subscript[A, 0]) and (Subscript[B, 0]) as its coefficients, and one represents far left Transmitted wave with coefficients (Subscript[A, N+1]) and (Subscript[B, N+1]), Then setting (Subscript[B, N+1]) = 0, and Solving for T = ((Subscript[A, N+1])/(Subscript[A, 0]))

(Subscript[A, N+1])==-I (β (Subscript[A, 0])+I (Subscript[A, 0])+β (Subscript[B, 0]) (E^(-2 I k (Subscript[x, 0]))))//TraditionalForm
(Subscript[B, N+1])==I β (Subscript[A, 0]) (E^(2 I k (Subscript[x, 0])))+I β (Subscript[B, 0])+(Subscript[B, 0])//TraditionalForm
0==I β (Subscript[A, 0]) (E^(2 I k (Subscript[x, 0])))+I β (Subscript[B, 0])+(Subscript[B, 0])//Simplify//TraditionalForm
AA=Solve[0==I β (Subscript[A, 0]) (E^(2 I k (Subscript[x, 0])))+I β (Subscript[B, 0])+(Subscript[B, 0]),(Subscript[A, 0])][[1]][[1]];
AA//TraditionalForm
T== ((Abs[((Subscript[A, N+1])/(Subscript[A, 0]))])^2)==ComplexExpand@Abs[((FullSimplify[-I (((β (Subscript[A, 0])+I (Subscript[A, 0])+β (Subscript[B, 0]) (E^(-2 I k (Subscript[x, 0])))))/(Subscript[A, 0]))/.(Subscript[A, 0])->(Subscript[A, 0])/.AA])^2)]//TraditionalForm

Now Recall M and lets check if T = (1/(|(Subscript[M, 22])((|)^2)))

"(Subscript[M, n])"==MatrixForm@Normal@CoefficientArrays[{(Subscript[A, n+1])/. H[n],(Subscript[B, n+1])/. H[n]},{(Subscript[A, n]),(Subscript[B, n])}][[2]]//TraditionalForm
"M =\!\(\*UnderoverscriptBox[\(\(\\ \)\(∏\)\), \(n = 0\), \(N\)]\) (Subscript[M, n])"//TraditionalForm
(Subscript[M, 22])==Part[M[n],2,2]//TraditionalForm
HoldForm@((1/((Abs[(Subscript[M, 22])])^2)) ) "(Overscript[=, ?])"  T//HoldForm//TraditionalForm
(1/((Abs[(Subscript[M, 22])])^2))==ComplexExpand@((Abs[(1/(Part[M[n],2,2]))])^2) == T//TraditionalForm

Bonus: Checking whether T + R = 1 holds

(Subscript[A, N+1])==-I (β (Subscript[A, 0])+I (Subscript[A, 0])+β (Subscript[B, 0]) (E^(-2 I k (Subscript[x, 0]))))//TraditionalForm
(Subscript[B, N+1])==I β (Subscript[A, 0]) (E^(2 I k (Subscript[x, 0])))+I β (Subscript[B, 0])+(Subscript[B, 0])//TraditionalForm
0==I β (Subscript[A, 0]) (E^(2 I k (Subscript[x, 0])))+I β (Subscript[B, 0])+(Subscript[B, 0])//Simplify//TraditionalForm
BB=Solve[0==I β (Subscript[A, 0]) (E^(2 I k (Subscript[x, 0])))+I β (Subscript[B, 0])+(Subscript[B, 0]),(Subscript[B, 0])][[1]][[1]];
BB//TraditionalForm
R==ComplexExpand@Abs[((FullSimplify[((-(((Subscript[A, 0]) β (E^(2 I k (Subscript[x, 0]))))/(β-I)))/(Subscript[A, 0]))/.(Subscript[B, 0])->(Subscript[B, 0])/.AA])^2)]//TraditionalForm

R+T==ComplexExpand@Abs[((FullSimplify[((-(((Subscript[A, 0]) β (E^(2 I k (Subscript[x, 0]))))/(β-I)))/(Subscript[A, 0]))/.(Subscript[B, 0])->(Subscript[B, 0])/.AA])^2)]+ComplexExpand@((Abs[(1/(Part[M[n],2,2]))])^2)==1//TraditionalForm

ploty[u_]:=Plot[ComplexExpand@((Abs[(1/(#))])^2)[[2]][[2]],{m,0,30},PlotRange->{{0,25},{0,1}},PlotLabel->n==u,ImageSize->Medium,PlotStyle->{{AbsoluteThickness[1.7],Black},{AbsoluteThickness[3.],Red}},MaxRecursion->5]&
" "
ns={0,1,5,10,20,100};
Ppl[u_]:=Plot[{ComplexExpand@((Abs[(1/(#))])^2)[[2]][[2]],1-ComplexExpand@((Abs[(1/(#))])^2)[[2]][[2]]},{m,0,10},PlotLegends->{"T","R"},PlotLabel->n==u,PlotRange->{{0,10},{0,1}},PlotStyle->{Black,Red},MaxRecursion->5]&

For[w={};i=0,i<6,i++,we=AppendTo[w,ploty[ns[[i+1]]]@wq[ns[[i+1]]]]]
exexe= Grid[{{w[[1]],w[[2]]},{w[[3]],w[[4]]},{w[[5]],w[[6]]}},Dividers->All]

Manipulate[Ppl[N]@wq[N],{N,0,5,1}]

Manipulate[ListLinePlot[Table[(1/(1+((((1/k)*Pi*ChebyshevU[n,(((-1))^n)]))^2))),{n,0,100,1}],PlotRange->All],{k,36,1200}]

ListLinePlot[Table[(1/(1+((((1/36)*Pi*ChebyshevU[n,(((-1))^n)]))^2))),{n,0,100,1}],PlotRange->All,PlotLabel->"T vs. N",PlotStyle->{Black}]

Button["HideInputs",NotebookFind[SelectedNotebook[],"Output",All,CellStyle];
FrontEndExecute[FrontEndToken["OpenCloseGroup"]]]