Using Mathematica To Find Reflection and Transmission Coefficients

In this notebook we'll demonstrate how to use Mathematica to compute reflection and transmission coefficients.  As an example, we'll compute and plot T and R for problem 3(a) of pset 9.

Using the continuity conditions for the finite well, one obtains the following equations that need to be solved simultaneously:
A (e^(-i k L/2))+B (e^(i k L/2))==C (e^(-i q L/2))+D (e^(i q L/2))
C (e^(i q L/2))+D (e^(-i q L/2))==F (e^(i k L/2))
k(A (e^(-i k L/2))-B (e^(i k L/2)))==q (C (e^(-i q L/2))-D (e^(i q L/2)) )
q (C (e^(i q L/2))-D (e^(-i q L/2)))==F k (e^(i k L/2)),
where k=((Sqrt[2 μ E])/ℏ)and q=((Sqrt[2 μ (E+V0)])/ℏ).  In setting up these equations I centered my well on x=0 (so the well goes from -L/2 to +L/2).  The final answers will not depend on where I place my well, but the initial equations may look a little different.

In solving this set of equations, Mathematica will actually be smart enough to realize that I have five unknowns (A, B, C, D, E) but only four equations.  It will therefore express everything in terms of ratios of coefficients.  However, it is easier to simply set A=1.  To solve the equations, we use the command Solve.  The basic syntax of Solve is Solve[{eqn1, eqn2, eqn3...}, {var1, var2, ...}], where eqn1 is the first equation (be sure to use double equal signs == in your equations), var1 is the first variable you want to solve for, and so on.

Solutions=Solve[{(E^(-I k L/2))+B (E^(I k L/2))==C (E^(-I q L/2))+D (E^(I q L/2)), C (E^(I q L/2))+D (E^(-I q L/2))==F (E^(I k L/2)),k((E^(-I k L/2))-B (E^(I k L/2)))==q (C (E^(-I q L/2))-D (E^(I q L/2)) ),q (C (E^(I q L/2))-D (E^(-I q L/2)))==F k (E^(I k L/2))}, {B, C, D, F}]

This gives us our solutions.  Next we need to do a little housekeeping.  Because of the way Mathematica works, we need to "tell" it that we want it to store certain results.  In what follows, the "/." command tells Mathematica to look at the vector called Solutions and to replace B with whatever B is in that vector.  (The [[1]] symbol tells Mathematica to look at the first element of the solutions vector, and is an annoying complication that arises because Mathematica decided to use double braces {{}} for the vector, so the entire vector is technically stored in the first element).

BoverA=B/.Solutions[[1]]

FoverA=F/.Solutions[[1]]

If we want, we can use the FullSimplify command to simplify the results (even though we really don't care about simplifying now that Mathematica's doing all the heavy algebraic lifting for us):

BoverA=FullSimplify[BoverA]

FoverA=FullSimplify[FoverA]

Getting the transmission and reflection coefficients can be somewhat yucky because Mathematica doesn't always know what quantities are real and what quantities are imaginary.  The follow commands seem to do the trick:

Transcoeff=FullSimplify[ComplexExpand[FoverA*Conjugate[FoverA]]]

Reflcoeff=FullSimplify[ComplexExpand[BoverA*Conjugate[BoverA]]]

We can now ask Mathematica to non-dimensionalize the expressions:

T[x_,g0_]=FullSimplify[Transcoeff/.k->((Sqrt[2 μ En])/ℏ)/.q->((Sqrt[2 μ (En+V0)])/ℏ)/.En->x V0/.L ->((g0 ℏ)/((Sqrt[2 μ])(Sqrt[V0])))];

R[x_,g0_]=FullSimplify[Reflcoeff/.k->((Sqrt[2 μ En])/ℏ)/.q->((Sqrt[2 μ (En+V0)])/ℏ)/.En->x V0/.L ->((g0 ℏ)/((Sqrt[2 μ])(Sqrt[V0])))];

Plotting the results for g0 = 20:

Geenought=20;
Plot[{T[x, Geenought], R[x,Geenought]}, {x, 0, 5}, PlotRange->{0, 1.1}]

Plotting the functions for general g0 using Mathematica's new dynamic capabilities:

Manipulate[Plot[{T[x, g0], R[x,g0]}, {x, 0, 5}, PlotRange->{0, 1.1}], {g0, 0, 50}]