Penetration through a Rectangular Barrier

Actually, some of the calculation is on the board, but this will do some messy stuff and make plots

Remove["Global`*"]

Write the expressions for the three wave functions

uI=cA Exp[I k x]+cB Exp[-I k x];
uII=cF Exp[q x]+cG Exp[-q x];
uIII=cC Exp[I k x];

Get the equations for the boundary conditions

eq1=(uI/.x->0)==(uII/.x->0)

eq2=(D[uI,x]/.x->0)==(D[uII,x]/.x->0)

eq3=(uII/.x->a)==(uIII/.x->a)

eq4=(D[uII,x]/.x->a)==(D[uIII,x]/.x->a)

Solve the first two equations for F and G

sol=Solve[{eq1,eq2},{cF,cG}]

Substitute these into the third and fourth equations and divide through by A

eq3x=((eq3/.sol/.{cB->rB cA,cC->rC cA})/.cA->1)[[1]]

eq4x=((eq4/.sol/.{cB->rB cA,cC->rC cA})/.cA->1)[[1]]

Solve for the ratios B/A and C/A

sol=Solve[{eq3x,eq4x},{rB,rC}];
BovrA=rB/.sol[[1]]
CovrA=rC/.sol[[1]]

Simplify[CovrA]

Write the transmission coefficient. This is from a hand calculation, using the result for C/A. It is tricky to take complex conjugates, for some reason.

T=1/(1+((k^2+q^2)/(2k q))^2Sinh[q a]^2);

Make a plot of T as a function of energy E, less than the barrier height V

kSub=k->Sqrt[2 m Enrg/ℏ^2];
qSub=q->Sqrt[2m (V-Enrg)/ℏ^2];

TEnrg=Simplify[T/.{kSub,qSub}]

TEpsilon=Simplify[TEnrg/.{Enrg->ϵ V}]

Vparam=16;
TPlot=Simplify[TEpsilon/.{V->Vparam ℏ^2/(2m a^2)},a>0]

Plot[TPlot,{ϵ,0,1},PlotStyle->{Thick,Black}]