THE AMMONIA MASER
				Maria Daghofer, Wolfgang von der Linden

Institut fuer Theoretische Physik
		Technische Universitaet Graz
		A-8010  Graz, Austria
		Projekt Multimediale Lehre des
		Zentralen Informatikdienstes der TU Graz

		2001-12-13

\(\(\(
\)\(\(\(Begin["ammonia`"]\) \)\(;\)\(\(Clear["ammonia`*"]\) \)\(;\)\(\ \)\)\)\)

Introduction: The (Subscript[NH, 3])molecule and its potential

The ammonia molecule ((Subscript[NH, 3])) consists of a pyramid with the N either above or below the plane of the H's:

Module[{pts, siz, lines, col, points, dl},
pts         =  {{0, 0, 0}, {1.2, 0, 
                0}, {0.5, (Sqrt[5/4]), 0}, {0.5, (Sqrt[5/4])/2, 2/(Sqrt[8])}};
siz         =  { .09, 0.09, 
              0.09,  .15};
col         =  {1, 1, 1,  .6};
points  = 
            Table[{Hue[col[[i]]], PointSize[siz[[i]]], 
                Point[pts[[i]]]}, {i, 1, 4}];
lines     =  
            Map[Line, 
              Flatten[Table[
                  Table[{pts[[i]], pts[[j]]}, {j, i + 1, 4}], {i, 1, 
                    3}], 1]];
dl             =  {Dashing[{0.05, 0.05}], Thickness[ .01], 
              Line[{{pts[[4, 1]], pts[[4, 2]], -pts[[4, 3]]}, 
                  pts[[4]]}]};

ppmol=Graphics3D[{PointSize[0.08], points, lines, dl},PlotRange -> {{- .1, 1.2}, {0, 
                      1.5}, {- .6,  .8}},Boxed -> False,Axes -> False,ImageSize  ->  300,ViewPoint->{1.786,  -2.497,  1.422}];
Show[ppmol];]

The three H-atoms can be above or below the nitrogen. Both configurations 
represent a minimum in the potential energy and the H-atoms can swing from 
one position to the other. This transition is used in an ammonia maser to 
produce microwaves. Alternatively, the H's can be taken as fixed and the 
nitrogen is then oscillating in a potential with two minima separated by a finite 
potential barrier.

Potential

The simplified potential which is going to be used , is computed and plotted:

Block[{g1, g2},
$TextStyle = {FontFamily -> "Times", 
              FontSize -> 16, FontWeight -> "Bold"};
g1 = 
            Graphics[{Thickness[ .02], 
                Line[{{-1, 100}, {-1, 0}, {- .25, 0}, {- .25, 
                      1}, { .25, 1}, { .25, 0}, {1, 0}, {1, 100}}]}, 
              PlotRange -> {{-1.5, 1.5}, {- .5, 2}}];
\(g2 = \ 
            Graphics[{Text[\*""(Subscript[V, 0]) →"", {\(-1.3\), 
                    1}], 
			Text["0 →", {\(-1.3\), 
                    0.05}], 
			Text["b", { .25, \(- .2\)}], 
			Text[" -b", {\(- .25\), \(- .2\)}], 
			Text["-a", {\(-1\), \(- .2\)}], 
			Text["a", {1, \(- .2\)}], 
			Text["← Δ →", {\(- .6\),  .2}], 
			Text["← Δ →", { .6,  .2}], 
			Text["↑ ∞", {1.2, 1.85}]}, \ \ \ PlotRange -> {{\(-1.5\), 
                    1.5}, {\(- .5\), 2.1}}]\);
pp = Show[{g1, g2}, ImageSize  ->  
                  240];]

Save plots to ps-file

The above plots can be saved to ps-files, you just have to uncomment ( i.e. delete (* and *) )the following cell. Replacing "EPS" by "GIF", "JPEG" ... allows you to save them to a different format You can of course set different directories.

Parameters

Parameters of the potential (see introduction):

Number of timesteps for the simulation:

Schroedinger equation

Parts of eigenfunctions

Solution for -a < x ≤ -b, using  ψ(-a)=0:

(Subscript[L, 1])[k_] = 
      Simplify[DSolve[{(Subscript[ψ, 1])''[x] + (k^2) (Subscript[ψ, 1])[x] == 
                0, (Subscript[ψ, 1])[-a] == 0}, (Subscript[ψ, 1])[x], x], 
          k ∈ Reals] /. {C[2] -> A , 
          C[1] -> B}

Solution for b ≤ x ≤ a, using  ψ(a)=0:

(Subscript[L, 3])[k_] = 
      Simplify[DSolve[{(Subscript[ψ, 3])''[x] + (k^2) (Subscript[ψ, 3])[x] == 
                0, (Subscript[ψ, 3])[a] == 0}, (Subscript[ψ, 3])[x], x], 
          k ∈ Reals] /. {C[2] -> C , 
          C[1] -> D}

Solution for -b ≤ x ≤ b:

(Subscript[L, 2])[k_] = 
      Simplify[DSolve[{(Subscript[ψ, 2])''[
                    x] + ((k^2) - (kv^2)) (Subscript[ψ, 2])[x] == 
                0}, (Subscript[ψ, 2])[x], x], {kv ∈ Reals, 
            k ∈ Reals}] /. {C[2] -> F, 
          C[1] -> G}

Eigenvalues and Eigenfunctions

As the potential is symmetric with respect to the y-axis, the  eigenfunctions must be either symmetric or antisymmetric. This can be  exploited to simplify the coefficients and gives

({{Symmetric eigenfunctions:, C  =  - A ,, G  =    F}, {Antisymmetric eigenfunctions:  , C  =    A ,, G  =  - F}})

Symmetric (subscript S) and antisymmetric (subscript A) eigenfunctions

\(\(ψ\_1[k_, x_] = \(ψ\_1[x] /. \(L\_1[k]\)[\([1]\)]\) /. 
          A -> A\&~/Sec[a\ k];\)
 (*D[ψ\_1[k, x], 
        x]*) \)
\(\(ψ\_\(3  S\)[k_, x_] = \(ψ\_3[x] /. \(L\_3[k]\)[\([1]\)]\) /. 
          C -> \(-A\&~\)/Sec[a\ k];\)\)
\(\(ψ\_\(3  A\)[k_, x_] = \(ψ\_3[x] /. \(L\_3[k]\)[\([1]\)]\) /. 
          C -> A\&~/Sec[a\ k];\)\)
\(\(ψ\_\(2  S\)[k_, x_] = 
        FullSimplify[\(\(ψ\_2[x] /. \(L\_2[k]\)[\([1]\)]\) /. 
              G -> F\) /. 
            F -> F\&~/2];\)
 (*D[ψ\_\(2  A\)[k, 
          x], x]*) \)
\(\(ψ\_\(2  A\)[k_, x_] = 
        FullSimplify[\(\(ψ\_2[x] /. \(L\_2[k]\)[\([1]\)]\) /. 
              G -> \(-F\)\) /. 
            F -> F\&~/2];\)\)
\(\(Print[
        TableForm[{{\*""(Subscript[ψ, 1])="", ψ\_1[k, 
                x], \*""(Subscript[ψ, 1])'="", 
              D[ψ\_1[k, x], 
                x]}, {\*""(Subscript[ψ, 2  S])="", ψ\_\(2  S\)[
                k, x], \*""(Subscript[ψ, 2  S])'="", 
              D[ψ\_\(2  S\)[k, x], 
                x]}, \ 
{\*""(Subscript[ψ, 2  A])="", ψ\_\(2  A\)[k, x], \*""(Subscript[ψ, 2  A])'="", 
              D[ψ\_\(2  A\)[k, x], 
                x]}, \ 
{\*""(Subscript[ψ, 3  S])="", ψ\_\(3  S\)[k, x], \*""(Subscript[ψ, 3  S])'="", 
              D[ψ\_\(3  S\)[k, x], 
                x]}, \ 
{\*""(Subscript[ψ, 3  A])="", ψ\_\(3  A\)[k, x], \*""(Subscript[ψ, 3  A])'="", 
              D[ψ\_\(3  A\)[k, x], x]}}]];\)\ \)

The wavefunction must be continuously differentiable, this  condition yields:

({{(Subscript[ψ, 1]) (-b)=(Subscript[ψ, 2  S])(-b) ,	(Subscript[(ψ'), 1]) (-b)=(Subscript[(ψ'), 2  S]) (-b)}, {
or:
}, {(Subscript[ψ, 1]) (-b)=(Subscript[ψ, 2  A])(-b) ,	(Subscript[(ψ'), 1]) (-b)=(Subscript[(ψ'), 2  A]) (-b).}})

Each pair of equations gives a system with the two independent  variables  \(\(\(A\)\(\ \)\)\&~\) and   \(F\&~\):

System for the symmetric eigenfunctions:

sysS  =  {(Subscript[ψ, 1])[k, -b]  -  (Subscript[ψ, 2  S])[
                k, -b]  ==  
            0, 
            (D[(Subscript[ψ, 1])[k, 
                      x], x]  /.  
                  x  ->  -b )  -   (D[(Subscript[ψ, 2  S])[k, 
                      x], x]  /.  x  ->  -b )  ==  
            0};
sysS // TableForm

System for the antisymmetric eigenfunctions:

sysA  =  {(Subscript[ψ, 1])[k, -b]  -  (Subscript[ψ, 2  A])[
                k, -b]  ==  
            0, 
            (D[(Subscript[ψ, 1])[k, 
                      x], x]  /.  
                  x  ->  -b )  -   (D[(Subscript[ψ, 2  A])[k, 
                      x], x]  /.  x  ->  -b )  ==  
            0};
sysA // TableForm

These systems of equations lead to the introduction of new, dimensionless  quantities:  kd  =  k (a - b),  β = b / (a - b)  and  kvd = (Subscript[k, V]) ·  (a - b).
 In order to give a non-trivial solution, the determinants of the  above systems must be 0.
  This gives us the following transcendental equations as eigenvalue equations:

Symmetric eigenfunctions:			tan(kd)  =  - \(TraditionalForm\`\_kd\/\@\(kvd\^2 - kd\^2\)\)coth[ β(Sqrt[(kvd^2) - (kd^2)])];
Antisymmetric eigenfunctions:		tan(kd)  = - \(TraditionalForm\`\_kd\/\@\(kvd\^2 - kd\^2\)\)tanh[ β(Sqrt[(kvd^2) - (kd^2)])]

This form of the determinants can be obtained from the systems of equations via some transformations:

Eigenvalue Equation for Symmetric Eigenfunctions

Coefficients of the symmetric system using the dimensionless quantities:

\(\(cmS[kd_]\  := \ 
        Simplify[{{Coefficient[sysS[\([1, 1]\)], A\&~], 
                Coefficient[sysS[\([1, 1]\)], \ 
                  F\&~]}, \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ {Coefficient[sysS[\([2, 1]\)], A\&~], 
                Coefficient[sysS[\([2, 1]\)], \ 
                  F\&~]}} /. {\((a - b)\) -> d, b -> β\ \ d, 
              k -> kd\ /\ d, kv -> kvd\ /\ d}, {k\  ∈ \ 
              Reals, \ k > 0}];\)\)
Print[
        ColumnForm[{"Coefficients of the symmetric system", "using the \
dimensionless quantities: "}], "    ", 
        MatrixForm[cmS[kd]]];
Print["The determinant of this system:  ", 
        Det[cmS[kd]]];
Print["This determinant is divided by the two coefficients stemming \
from the second equation and multiplied by kd:"];
symmdet[kd_] := 
        Simplify[
          kd Det[cmS[
                  kd]]/(d*cmS[kd][[2, 1]]*cmS[kd][[2, 
                      2]]), {d ϵ Reals, 
            d > 0}];
symmdet[kd]

Eigenvalue Equation for Antisymmetric Eigenfunctions

\(\(cmA[kd_]\  := \ 
        Simplify[{{Coefficient[sysA[\([1, 1]\)], A\&~], 
                Coefficient[sysA[\([1, 1]\)], \ 
                  F\&~]}, \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ {Coefficient[sysA[\([2, 1]\)], A\&~], 
                Coefficient[sysA[\([2, 1]\)], \ 
                  F\&~]}} /. {\((a - b)\) -> d, b -> β\ \ d, 
              k -> kd\ /\ d, kv -> kvd\ /\ d}, {k\  ∈ \ 
              Reals, \ k > 0}];\)\)
Print[
        ColumnForm[{"Coefficients of the antisymmetric system", "using \
the dimensionless quantities: "}], "    ", 
        MatrixForm[cmA[kd]]];
Print["The determinant of this system:  ", 
        Det[cmA[kd]]];
Print["This determinant is divided by the two coefficients stemming \
from the second equation and multiplied by kd:"];
asymmdet[kd_] := 
        Simplify[
          kd Det[cmA[
                  kd]]/(d*cmA[kd][[2, 1]]*cmA[kd][[2, 
                      2]]), {d ϵ Reals, 
            d > 0}];
asymmdet[kd]

Determination of the Eigenvalues

The zeroes of the functions symmdet and asymmdet give the momenta of the eigenfunctions, the eigenenergies are then obtained via k =  ((Sqrt[2 m E])/ℏ). Approximate values can be found by a graphical method.
For given values of  V0, a  and b  (thus for a given value of kvd), the two summands of symmdet are plotted versus kd. 
Intersections of the curves give the zeroes kd.  The same applies to faa.
Here, numerical values are obtained by FindRoot[]; these are stored in the lists  ks and ka.

kvd := (Sqrt[2  V0])(a - b);β := b/(a - b);
\(Block[{a, b, V0, nn, ϵ, k0, 
        k1, n, km, fs1, fs2, fa2}, 
a = a0; b = b0; 
      V0 = V00; 
ϵ\ \ \  = \  .00000000001; 
nn\  = \ 
        Floor[\(2\ kvd\  + \ π\)\/\(2\ π\)]; 
fs1 = \(symmdet[kd]\)[\([2]\)]; fs2 = \(-\(symmdet[kd]\)[\([1]\)]\); 
      fa2 = \(-\(asymmdet[kd]\)[\([2]\)]\); 
Plot[{fs1, 
          fs2, fa2}, {kd, 0, kvd}, PlotRange -> {{0, kvd}, {\(-5\), 5}}, 
        AxesLabel -> {"kd", ""}, 
        PlotStyle -> {Hue[ .6], Hue[1], Hue[0.3]}, 
        TextStyle -> {FontWeight -> Plain}, 
        ImageSize -> 300]; 
ks\  = 
        Flatten[\ 
          Table[\ k1\  = \ Min[n\ π, kvd] - ϵ; \ 
            k0\  = \ n\ π - π/2 + ϵ; \ 
            km\  = \ \((k0 + k1)\)/2; 
			If[symmdet[k0]*symmdet[k1] < 0, 
              kd /. FindRoot[
                  symmdet[kd] == 0, {kd, km, k0, 
                    k1}], {}]
\ \ \ \ \ \ \ \ \ \ , {n, 1, 
              nn}]]; 
ka\  = \ 
        Flatten[Table[\ k1\  = \ Min[n\ π, kvd] - ϵ; \ 
            k0\  = \ n\ π - π/2 + ϵ; \ 
            km\  = \ \((k0 + k1)\)/2; 
			If[
              asymmdet[k0]*asymmdet[k1] < 0, 
              kd /. FindRoot[
                  asymmdet[kd] == 0, {kd, km, k0, 
                    k1}], {}]
, {n, 1, 
              nn}]];
]\)

The intersections of the blue curves for the tangent with the red (belonging to symmetric wavefunctions) 
and green  (belonging to antisymmetric wavefunctions)  curves give the possible values of kd for momenta of the eigenfunctions.
Their numerical values and the eigenenergies are given in the lists below:
Symmetric solutions:

Table[{"k:", ks[[ik]], "    E:", ks[[ik]]^2/2}, {ik, 
          1, Length[ks]}] // TableForm

Antisymmetric solutions:

Table[{"k:", ka[[ik]], "    E:", ka[[ik]]^2/2}, {ik, 
          1, Length[ka]}] // TableForm

Plot of eigenfunctions

For the k values obtained , we can choose \(TraditionalForm\`\(\(A\)\(\ \)\)\&~\) and \(TraditionalForm\`F\&~\) to define continuos  wave-functions , the constants are finally determined by the normalization.

\(φs[i_, x_] := \ \ Which[x < \(-1\) - β, 0, \ 
        x < \ \(-β\), \ 
        As\&~[\([i]\)]\ \ Sin[ks[\([i]\)] \((1 + β + x)\)], \ 
        x\  < \ β, 
        Fs\&~[\([i]\)]\ Cosh[\@\(\(kvd\^2\)\(-\)\(\(\(ks[\([i]\)]\)\(\ \)\)\^2\)\(\ \)\)\ x\ ], 
			\ \ \ \ \ \ \ \ \ \ \ \ \ x < 1 + β, As\&~[\([i]\)]\ Sin[ks[\([i]\)] \((1 + β - x)\)], \ 
        x ≥ \ 1 + β, 0]\)
\(φa[i_, x_] := \ 
      Which[x < \(-1\) - β, 0, \ x < \ \(-β\), \ 
        Aa\&~[\([i]\)]\ \ Sin[ka[\([i]\)] \((1 + β + x)\)], 
        x\  < \ β, 
        Fa\&~[\([i]\)]\ Sinh[\@\(\(kvd\^2\)\(-\)\(\(\(ka[\([i]\)]\)\(\ \)\)\^2\)\(\ \)\)\ x], 
				x < 
          1 + β, \(-Aa\&~[\([i]\)]\)\ Sin[
            ka[\([i]\)] \((1 + β - x)\)], \ 
        x ≥ \ 1 + β, 0]\)

\(Block[{b, a, V0, i, x}, 
b = b0; a = a0; 
      V0 = V00; 
As\&~ = 
        Cosh[\@\(\(kvd\^2\)\(-\)\(ks\^2\)\(\ \)\)\ β]; 
      Fs\&~ = \ Sin[ks]; 
Zs = 
        Table[NIntegrate[Abs[φs[i, x]]\^2, {x, \(-a\), a}], {i, 1, 
            Length[ks]}]; 
As\&~ = 
        Cosh[\@\(\(kvd\^2\)\(-\)\(ks\^2\)\(\ \)\)\ β]/\@Zs; 
      Fs\&~ = \ Sin[ks]/\@Zs; 
Aa\&~ = 
        Sinh[\@\(\(kvd\^2\)\(-\)\(ka\^2\)\(\ \)\)\ β]; 
      Fa\&~ = \ \(-Sin[ka]\); 
Za = 
        Table[NIntegrate[Abs[φa[i, x]]\^2, {x, \(-a\), a}], {i, 1, 
            Length[ka]}]; 
Aa\&~ = 
        Sinh[\@\(\(kvd\^2\)\(-\)\(ka\^2\)\(\ \)\)\ β]/\@Za; 
      Fa\&~ = \ \(-Sin[ka]\)/\@Za;
]\)

\(Block[{b, a, V0, ik, x}, 
b = b0; a = a0; 
      V0 = V00; 
ϕls\  = \ 
        Table[φs[ik, x], {ik, 
            Length[ks]}]; 
ϕla\  = \ 
        Table[φa[ik, x], {ik, 
            Length[ka]}]; 
pll\  = \ 
        Table[Hue[ik/Length[ks]], {ik, 0, 
            Length[ks]\  - 1}]; 
psymm = 
        Plot[Evaluate[ϕls], {x, \(-a\), a}, 
          AxesLabel -> {\*""\!\(x\/\(a - b\)\)"", ""}, 
          PlotLabel -> "symmetric", PlotStyle -> pll, 
          TextStyle -> {FontWeight -> 
                Plain}]; 
pasymm = 
        Plot[Evaluate[ϕla], {x, \(-a\), a}, 
          AxesLabel -> {\*""\!\(x\/\(a - b\)\)"", ""}, 
          PlotLabel -> "antisymmetric", PlotStyle -> pll, 
          TextStyle -> {FontWeight -> 
                Plain}];
]\)

Dynamics

\(
\(dyn[ca_, cs_, tmax_, Nt_] := 
        Block[{x, na = Length[ca], ns = Length[cs], j, p1, it, 
            dt = tmax/Nt}, 
F[x_, t_] := 
            Sum[Exp[I\ t\ \((ks[\([j]\)])\)\^2/
                        2]\ cs[\([j]\)]*φs[j, x], {j, 1, ns}] + 
              Sum[Exp[I\ t\ \((ka[\([j]\)])\)\^2/2] 
                  ca[\([j]\)]*φa[j, x], {j, 1, 
                  na}]; 
p1 = 
            Plot[Abs[F[x, 0]]\^2, {x, \(-a\), a}, 
              PlotRange -> {{\(-a\), a}, {0, 2}}, 
              DisplayFunction -> Identity]; 
ptable = 
            Table[Show[{p1, 
Plot[
                    Abs[F[x, dt\ it]]\^2, {x, \(-a\), a}, 
                    PlotStyle -> {Hue[1]}, 
                    PlotRange -> {{\(-a\), a}, {0, 2}}, 
                    DisplayFunction -> Identity]}, 
                AxesLabel -> {"x", "ρ(x)"}, 
                TextStyle -> {FontWeight -> Plain}, 
                DisplayFunction -> $DisplayFunction], {it, 1, 
                Nt}]; 
	
          SelectionMove[EvaluationNotebook[], After, 
            EvaluationCell]; 
	
          SelectionMove[EvaluationNotebook[], Next, 
            CellGroup]; 
	
          FrontEndTokenExecute["OpenCloseGroup"];
];\)\)
start = 
        ButtonBox[
          StyleBox["Run Animation", "Title", 
            FontFamily -> "Helvetica"], 
          ButtonFunction :> 
            FrontEndExecute[{SelectionMove[EvaluationNotebook[], Before, 
                    ButtonCell]
SelectionMove[
                    EvaluationNotebook[], Previous, 
                    CellGroup]
SelectionAnimate[
                    SelectedNotebook[], 
                    AnimationDisplayTime -> simtime]}], 
          ButtonEvaluator -> Automatic, 
          Active -> True];

This simulation shows the dynamics of a superposition of the first two 
eigenstates (i.e. of one symmetric and one antisymmetric eigenstate):

Block[{a = a0, b = b0, V0 = V00, ca, 
          cs}, 
  
ca = {1/N[(Sqrt[2])], 
            0}; 
cs = {1/N[(Sqrt[2])], 0}; 
dyn[
          ca, cs, 4 π/((((ka[[1]]))^2) - (((ks[[1]]))^2)), 
          Nt0]];
\(start /. simtime -> \ 0.1 // DisplayForm	\)

If you want to save this animation as an animated GIF file, you have to 
uncomment the following cell. Just delete the (* and the *), saving may take 
fairly long.

The period of this inter-well oscillation is (Subscript[τ, 1])=\(TraditionalForm\`\(4  π\)\/\(k\_\(1, S\)\%2 - k\_\(1, A\)\%2\)\):

\(\(\(
\)\(Print[\*""\!\(τ\_1\)= "", 
        4\ π/\((\((ka[\([1]\)])\)\^2 - \((ks[\([1]\)])\)\^2)\)];\)\)\)

This simulation shows the dynamics of a superposition of the first two 
symmetric (i. e. the first and the third) eigenstates:

\(Block[{a = a0, b = b0, V0 = V00, ca, 
        cs}, 
ca = {0, 
          0}; 
cs = {1/N[\@2], 
          1/N[\@2]}; 
dyn[ca, cs, 
        4\ π/\((\((ks[\([2]\)])\)\^2 - \((ks[\([1]\)])\)\^2)\), 
        Nt0]]; \(\(start /. 
        simtime -> \ 0.1\)\(//\)\(DisplayForm\)\(	\)\)\)

Nt2 = 
        Round[2 ((((ks[[2]]))^2) - (((ks[[1]]))^2))/((((ka[[1]]))^2) - (((ks[[1]]))^2))];
\(\(Print[\*""(Subscript[τ, 2])= "", 
        4\ π/\((\((ks[\([2]\)])\)\^2 - \((ks[\([1]\)])\)\^2)\), \*""        (Subscript[τ, 1])/(Subscript[τ, 2])= "", \((\((ks[\([2]\)])\)\^2 - \((ks[\([1]\)])\)\^2)\)/\((\((ka[\([1]\)])\)\^2 - \((ks[\([1]\)])\)\^2)\)];\)\)

The timescale of these two motions can be seen in this simulation, which  shows a superposition of the first three eigenstates:

● Caution: The number of timesteps is roughly 2(Subscript[τ, 1])/(Subscript[τ, 2]) which can be rather large!

\(Block[{a = a0, b = b0, V0 = V00, ca, 
        cs}, 
cs = {0.5, 
          0.5}; 
ca = {1/\@2, 0}; 
dyn[
        ca, cs, 4\ π/\((\((ka[\([1]\)])\)\^2 - \((ks[\([1]\)])\)\^2)\), 
        Nt2]]; \(\(start /. 
        simtime -> \ 0.1\)\(//\)\(DisplayForm\)\(	\)\)\)

Exercises

1.  Higher or Lower Steps

You may choose a higher value for V00, which corresponds to a higher step in the center separating the two potential wells.
If this value is high enough you get 6 in place of 4 eigenvalues.  If you take V00 sufficiently small you get only 
2 in place of 4 eigenvalues. To do this you change the value of  V00 in the field for the input parameters and evaluate the
notebook again.

2.  Finding Critical Values of Step Height

By using the method suggested in 1. you may find the critical values of  V00 separating the domains giving 2, 4, 6,...
eigenvalues.