(*^
::[	Information =

	"This is a Mathematica Notebook file.  It contains ASCII text, and can be
	transferred by email, ftp, or other text-file transfer utility.  It should
	be read or edited using a copy of Mathematica or MathReader.  If you 
	received this as email, use your mail application or copy/paste to save 
	everything from the line containing (*^ down to the line containing ^*)
	into a plain text file.  On some systems you may have to give the file a 
	name ending with ".ma" to allow Mathematica to recognize it as a Notebook.
	The line below identifies what version of Mathematica created this file,
	but it can be opened using any other version as well.";

	FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2";

	MacintoshStandardFontEncoding; 
	
	fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8,  24, "Times"; 
	fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6,  18, "Times"; 
	fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6,  14, "Times"; 
	fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20,  18, "Times"; 
	fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15,  14, "Times"; 
	fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12,  12, "Times"; 
	fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  14, "Times"; 
	fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, "Times"; 
	fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-5,  12, "Courier"; 
	fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5,  12, "Courier"; 
	fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5,  12, "Courier"; 
	fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5,  12, "Courier"; 
	fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5,  12, "Courier"; 
	fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, groupLikeGraphics, M7, l34, w405, h197,  12, "Courier"; 
	fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic,  10, "Geneva"; 
	fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = leftheader, inactive, L2,  12, "Times"; 
	fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7,  12, "Times"; 
	fontset = leftfooter, inactive, L2,  12, "Times"; 
	fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, "Times"; 
	fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	paletteColors = 128; currentKernel; 
]
:[font = section; inactive; Cclosed; noPageBreak; preserveAspect; cellOutline; backColorRed = 58981; backColorGreen = 58981; backColorBlue = 58981; startGroup]
Setup 
:[font = input; noPageBreak; preserveAspect]
Clear["Global`*"];
:[font = subsection; inactive; noPageBreak; preserveAspect; cellOutline; startGroup]
4. Complex Conjugate Rule 
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
We will find it convenient to introduce our own definition for \com{congujate} rather than use the built in \com{Conjugate}. 
:[font = input; noPageBreak; preserveAspect]
conjugate::usage = 
" A simple method of computing the conjugate of an object
 which is explicitly Complex";

conjugateRule     = {Complex[re_,im_]:>Complex[re,-im]};
conjugate[exp__] := exp /. conjugateRule;
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
Note the use of the function ensures we do not apply this in an endless loop.
:[font = input; noPageBreak; preserveAspect; endGroup; endGroup]
Protect[{conjugate ,conjugateRule}];
:[font = section; inactive; Cclosed; noPageBreak; preserveAspect; cellOutline; backColorRed = 58981; backColorGreen = 58981; backColorBlue = 58981; startGroup]
Example 1: Particle Propagating Towards a Rectangular Potential  
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; cellOutline]
Suppose a particle with energy En= Energy[n]>0 propagates from the left (x < -a)  towards a rectangular well or  barrier described by: V(x) =  V0   if -a<x<a, and  V(x) =0   otherwise.  
 
If V0 > 0, the potential describes a barrier, and if V0 < 0 the potential describes a well. The wave to the left (x < -a) of the potential has an incoming and reflected component. We set the amplitude of the incoming wave to one and let R be the amplitude of the reflected wave. The wave to the right (x > a) of the potential has only a transmitted wave; let its amplitude be T. The wave in the potential can be described by two amplitudes, call them A and B.  

a) Write the general form of the wave function in the  left, center, and right of the well.   
b) Use the boundary condition to evaluate the wave amplitudes. 
c) The coefficients of  transmission (T) and reflection (R) are defined as the ratio of the transmitted and reflected currents to the incident current.  Find the transmission and reflection coefficients.  Verify that T+R=1.  
d) Plot R and T vs. En  for the case of a potential well (V0<0)  and for the case of a potential barrier (V0>0).   
e) Use Plot3D to plot the reflection and transmission curves as a function of the energy and potential.  
f) Plot the real part of the wave for the complete range of x.
g) Animate the time dependence of the real part of the wave as it propagates through  a barrier  (V0>0, 0<En<V0).          
:[font = subsubsection; inactive; noPageBreak; preserveAspect; startGroup]
Remarks and Outline:
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; endGroup]
The general solution of Schrodinger's equation in a constant potential is given by a plane wave.  The wave solutions to the left, right, and in the center of the potential must be treated separately.  The  boundary conditions determine the wave function amplitudes.   The coefficients of  transmission (TT) and reflection (RR) follow  from computing the  flux  of the wave function.  Plots of these coefficients  reveal  the resonant  behavior. The reader is encouraged to vary the energy, potential, and well size. We define a graphic function to join the wave function in the three distinct potential regions. Finally the animation follows from forming a table of wave plots for different values of the time.
:[font = subsubsection; inactive; noPageBreak; preserveAspect; startGroup]
Required Packages
:[font = input; noPageBreak; preserveAspect; endGroup]
Needs["Algebra`Trigonometry`"]
Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`ParametricPlot3D`"]
:[font = subsubsection; inactive; noPageBreak; preserveAspect; startGroup]
Solution
:[font = input; noPageBreak; preserveAspect]
Clear["Global`*"];
:[font = input; noPageBreak; preserveAspect]
Off[General::spell ];
Off[General::spell1];
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; cellOutline]
a) Write the general form of the wave function in the  left, center, and right of the well.   
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; startGroup]
The general solution for Schrodinger's equation with a constant potential is given bythe following. For x < -a  there is an incident wave traveling to the right, and a reflected wave traveling to the left.  The amplitude for the incident wave is set equal to one without  loss of generality, and the amplitude for the reflected wave is called R. The form of the wave is
:[font = input; noPageBreak; preserveAspect]
psiL[x_]= E^(I k1 x) + E^(-I k1 x) R
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
For x > +a there is only a transmitted wave with amplitude T. 
:[font = input; noPageBreak; preserveAspect]
psiR[x_]= E^(I k1 x) T
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
The wave function in the central region  (-a<x<a) is 
:[font = input; noPageBreak; preserveAspect]
psiW[x_]= B E^(-I k2 x) + A E^(I k2 x)
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
The expressions for k1 and k2 follow from kRule  
:[font = input; noPageBreak; preserveAspect]
kRule= { k1->Sqrt[2 m (En  )]/hbar
        ,k2->Sqrt[2 m (En-V)]/hbar
       } //PowerExpand
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; endGroup]
where the wave number is implictly defined by  k^2 hbar^2 /(2m) = (En-V).
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; cellOutline]
b) Use the boundary condition to evaluate the wave amplitudes. 
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; startGroup]
The boundary conditions require the  wave function and its derivative  be continuous across the boundaries. Applying these conditions at  x=a and x=-a  yields  four equations  
:[font = input; noPageBreak; preserveAspect]
eq1=  {    (psiL[ x] - psiW[x]==0)   /. {x->-a}, 
           (psiW[ x] - psiR[x]==0)   /. {x->+a},
           (psiL'[x] -psiW'[x]==0)   /. {x->-a},
           (psiW'[x] -psiR'[x]==0)   /. {x->+a} }
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; endGroup]
We solve these equations for the amplitudes {A,B,R,T}. 
(This calculation takes a few minutes; it may be necessary to break the solution up into smaller steps on your particular machine.)
;[s]
3:0,0;56,1;187,0;188,-1;
2:2,16,12,Times,0,14,0,0,0;1,13,9,Times,2,12,0,0,0;
:[font = input; noPageBreak; preserveAspect]
ABRTrule=( Solve[eq1,{A,B,R,T}][[1]] 
               //Map[Together,#,{2}]&
               //ComplexToTrig 
               //Simplify
             )
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; cellOutline]
c) The coefficients of  transmission (T) and reflection (R) are defined as the ratio of the transmitted and reflected currents to the incident current.  Find the transmission and reflection coefficients.  Verify that T+R=1.  
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
The coefficients of  transmission and reflection are defined as the ratio of the transmitted and reflected currents as compared to the incident current. The vector currents are in the x direction so we only consider the first component. The currents follow from applying the user-defined function flux to the wave functions.  The x-component of the incident, reflected, and transmitted  currents are:
:[font = input; noPageBreak; preserveAspect; startGroup]
SetCoordinates[Cartesian[x,y,z]];   (* Used by Grad *)

flux[psi_]:= 
 ( hbar/(2 I m) 
   ( (psi//conjugate) Grad[psi] - Grad[psi//conjugate] (psi) ) 
 )
:[font = input; noPageBreak; preserveAspect; endGroup]
incFlux= flux[  psiL[x][[1]] ] 
:[font = input; noPageBreak; preserveAspect]
refFlux= flux[  psiL[x][[2]] ] 
:[font = input; noPageBreak; preserveAspect]
tranFlux= flux[  psiR[x] ] 
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
We find the reflection and transmission coefficients to be:
:[font = input; noPageBreak; preserveAspect]
RR=  refFlux[[1]]/incFlux[[1]] //Abs
TT= tranFlux[[1]]/incFlux[[1]] //Abs
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
Note, {RR,TT} are defined to be positive definite. We see that with the normalization we have chosen for the incident wave, that RR=R R*, and TT=T  T^*. In terms of the initial variables, we find:
:[font = input; noPageBreak; preserveAspect]
{TT,RR}= (( {T Conjugate[T], R Conjugate[R]}
         /.ABRTrule 
         /.{Conjugate->conjugate} 
         //Together 
         //Simplify
       ) //Map[Collect[#,{Cos[4 a k2]}]&,#,{3}]&
         //Simplify
       )
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
Notice we verify probability conservation.
:[font = input; noPageBreak; preserveAspect]
TT+RR==1 //ExpandAll//Simplify
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; cellOutline]
d) Plot R and T vs. En  for the case of a potential well (V0<0)  and for the case of a potential barrier (V0>0). 
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; startGroup]
We  plot RR and TT as a function of En for the parameters   
:[font = input; noPageBreak; preserveAspect]
values={hbar->1,m->1,a->1};
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
Consider the plots for the case that V<0  (potential well) and for the case V>0 (potential barrier). The graphics for  {RR,TT follow from the  function,  
:[font = input; noPageBreak; preserveAspect]
Clear[plotOperator];
plotOperator[Potential_, EnRange_] := 
Plot[ {RR,TT} /.kRule /.V->Potential /.values //Evaluate
  ,EnRange 
  ,AxesLabel->{"En","RR, TT"}
  ,PlotStyle->{{RGBColor[1,0,0],Dashing[{1,0.0}]}
              ,{RGBColor[0,0,1],Dashing[{0.02,0.02}],Thickness[0.008]}
              }
];  
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; endGroup]
The function plotOperator generates the graphics for {RR,TT}  as a function of energy and potential.  The RR coefficient is denoted by a  solid curve, and the TT coefficient is denoted by a dashed curve.  
We first consider the case of a potential well with V=-40 and plot the coefficients for En between 0 and 40.
:[font = input; noPageBreak; preserveAspect]
plotOperator[-40,{En,0,40}];
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
The dashed curve is the TT and the solid curve is RR. The striking decrease in the reflection coefficient with increasing energy was noted for gases by  Ramsauer in 1920,  and independently by Townsend and Bailey later.

Consider a potential barrier with V=1, and let the values of En range from 0.01 to 0.9.
:[font = input; noPageBreak; preserveAspect]
plotOperator[1,{En,0.01,0.99}];
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; cellOutline]
e) Use Plot3D to plot the reflection and transmission curves as a function of the energy and potential.  
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; startGroup]
We use Plot3D to illustrate the reflection and transmission curves for a range of potential and energy values. Applying Plot3D to RR and TT and plotting the curves as a function of energy and potential, we get  
:[font = input; noPageBreak; preserveAspect]
pt1a= 
Plot3D[TT /.kRule /.ABRTrule /.values /.a->1 //Evaluate
      ,{En,.01,10}
      ,{V,-20,-.1}
      ];
:[font = input; noPageBreak; preserveAspect]
pt1b= 
Plot3D[RR /.kRule /.ABRTrule /.values /.a->1 //Evaluate
      ,{En,.01,10}
      ,{V,-20,-.1}
      ];
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; endGroup]
Resonance effects appear as ripples in the diagram. 
 
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; cellOutline]
f) Plot the real part of the wave for the complete range of x.
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; startGroup]
We plot the real part of the wave for the complete range of x. We choose a=2, m=1 and hbar=1, and define plotOp
:[font = input; noPageBreak; preserveAspect]
plotOp[wave_,Energy_,V0_,thickness_,xrange_]:= 
  Plot[  wave /.ABRTrule 
              /.kRule
              //.{m->1,hbar->1,En->Energy,a->2,V->V0} 
              //Re 
              //Evaluate
    ,xrange
    ,PlotStyle-> Thickness[thickness] 
    ,AxesOrigin->{0,0} 
    ,DisplayFunction->Identity
];
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
The function plotOp generates the graphics for the wave with energy Energy and potential V0. The variable thickness allows us to distinguish the  wave to the  left, center, and right of the barrier. We combine graphics  with the function
:[font = input; noPageBreak; preserveAspect]
Clear[wavePlot1];
wavePlot1[En_,V0_,xmin_:-12,xmax_:12] :=
  Show[ {Graphics[{Hue[0.20,.50,1],Rectangle[{-2,-10^9},{+2,+10^9}]}] (* do this first *)
        ,plotOp[psiL[x],En,V0,0.008,{x,xmin,-2}]
        ,plotOp[psiR[x],En,V0,0.008,{x,2,xmax} ]
        ,plotOp[psiW[x],En,V0,0.004,{x,-2,2}   ]
        }
    ,DisplayFunction-> $DisplayFunction ];
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
The outline of the well is given in the command Graphics.
Displaying the real part of the wave for V0=-10 and En=1, we obtain 
:[font = input; noPageBreak; preserveAspect; startGroup]
wavePlot1[1,-10];
:[font = section; inactive; noPageBreak; preserveAspect; fontColorRed = 65535; endGroup; endGroup]
Your Challenge: Experiment with other combinations:
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; cellOutline]
g) Animate the time dependence of the real part of the wave as it propagates through  a barrier  (V0>0, 0<En<V0).          
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; startGroup]
Let us animate the real and imaginary part of the wave as it propagates through  the barrier (0<En<V0). To get the time dependence for the wave, we multiply the wave by 
Exp[-i E_n t/hbar] with hbar=1.  The time dependence of a wave with energy Energy  is described by the function,
:[font = input; noPageBreak; preserveAspect]
Clear[timePlot];
timePlot[time_,wave_,Energy_,V0_,xrange_]:= 
  Plot[{Exp[-I Energy time] wave
             /.ABRTrule
             /.kRule
             //.{m->1,hbar->1,En->Energy,a->2,V->V0}
             //Re  
       ,Exp[-I Energy time] wave
             /.ABRTrule
             /.kRule
             //.{m->1,hbar->1,En->Energy,a->2,V->V0}
             //Im  
        }//Evaluate  
    ,xrange
    ,Axes->None
    ,PlotStyle->{{RGBColor[1,0,0],Thickness[0.006]}
                ,{RGBColor[0,0,1],Thickness[0.006]}
                }
    ,DisplayFunction->Identity   
 ]; 
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
Let's just test this function:
:[font = input; noPageBreak; preserveAspect]
timePlot[1.234,psiL[x],1.0,1.1,{x,-10,-2}] //Show[#,DisplayFunction->$DisplayFunction]&;
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
Let us make a function (showit) to animate a sequence of plots given and Energy an V0.
:[font = input; noPageBreak; preserveAspect]
Clear[showit];
showit[EnergyIn_,V0In_,framesIn_:1]:= 
Module[{Energy,V0,tnumber,timestep},
 Energy=EnergyIn;
 V0=V0In;
 tnumber=framesIn;
 timestep=2 Pi/(tnumber);
 Do[ Show[ {Graphics[{Hue[0.20,.50,1],Rectangle[{-2,-10^9},{+2,+10^9}]}] (* do this first *)
           ,timePlot[time,psiL[x],Energy,V0,{x,-10,-2}]
           ,timePlot[time,psiW[x],Energy,V0,{x,-2,  2}]
           ,timePlot[time,psiR[x],Energy,V0,{x, 2, 10}]
           }
        ,DisplayFunction->$DisplayFunction   
        ,PlotRange->{{-10,10},{-2,2}}
        ] 
   ,{time,0,2 Pi-timestep,timestep} ]  

]
           
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
Here is a sample:
:[font = input; noPageBreak; preserveAspect]
showit[1.0,1.1]
:[font = section; inactive; noPageBreak; preserveAspect; fontColorRed = 65535]
Your Challenge: Experiment with other combinations:
:[font = input; noPageBreak; preserveAspect]
showit[1.0,-1.1]
:[font = input; noPageBreak; preserveAspect]
showit[1.0,+10.]
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; endGroup; endGroup; endGroup]
The number of animation frames is determined by  tnumber. The reader is encouraged to choose different values for the energy and potential.  Compare the results with the plot of TT and RR vs. energy. 
:[font = section; inactive; Cclosed; noPageBreak; preserveAspect; cellOutline; backColorRed = 58981; backColorGreen = 58981; backColorBlue = 58981; startGroup]
Example 2: Particle bound in an Infinite Potential Well  
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
Suppose a particle is bound by the infinite potential well  V(x) = If[0 < x < a; 0,Infinity].
        
a) Solve Schrodinger's equation for the general wave function when V=0.  
Use the boundary conditions at the walls to determine the eigenfunctions.  Normalize the eigenfunctions. 

b) Explicitly verify that the first four wave functions are  orthonormal, and graph  the wave functions.   

c) Determine the energy eigenvalues. 
:[font = subsubsection; inactive; noPageBreak; preserveAspect; startGroup]
Remarks and Outline:
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; endGroup]
This simple one-dimensional boundary problem illustrates the eigen solutions of Schrodinger's equation.  The  boundary conditions require the wave function to vanish at x=0 and a. 
:[font = subsubsection; inactive; noPageBreak; preserveAspect; startGroup]
Required Packages
:[font = input; noPageBreak; preserveAspect; endGroup]
Needs["Algebra`Trigonometry`"]
Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`ParametricPlot3D`"]
:[font = subsubsection; inactive; noPageBreak; preserveAspect; startGroup]
Solution: 
:[font = input; noPageBreak; preserveAspect]
Clear["Global`*"];
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; cellOutline]
a) Solve Schrodinger's equation for the general wave function when V=0.  
Use the boundary conditions at the walls to determine the eigenfunctions.  Normalize the eigenfunctions. 
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; startGroup]
Schrodinger's equation is given by:
:[font = input; noPageBreak; preserveAspect]
eq1= -En wave[x] -  wave''[x] hbar^2/(2 m) ==0
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
Here, we have taken the potential to zero. 
The solution follows from  DSolve:
:[font = input; noPageBreak; preserveAspect]
dsol= DSolve[eq1,wave,x][[1]]
          
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
We find it convenient to define the wave number, k:
:[font = input; noPageBreak; preserveAspect]
eq2= {k == Sqrt[2 En m]/hbar};

eSol= Solve[ eq2,En][[1]]
:[font = input; noPageBreak; preserveAspect]
kSol= Solve[ eq2,k][[1]]
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
In terms of k, the solution is:
:[font = input; noPageBreak; preserveAspect; endGroup]
dsol= dsol /.eSol //PowerExpand;
wave[x] /.dsol //Simplify
wave[x] /.dsol //ComplexToTrig //Simplify
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; cellOutline]
b) Explicitly verify that the first four wave functions are  orthonormal, and graph  the wave functions.   
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; startGroup]
The constants  {C[1], C[2]} are determined by the boundary conditions at the  infinite barrier, and by normalizing the wave function.  Requiring the wave function  vanish at  x={0,a}  dictates that either:
           1) {C[1], C[2]}={0,0}, a trivial solution, or 
           2) Sin[a k] = 0, a more interesting solution. 
:[font = input; noPageBreak; preserveAspect]
eq3={ wave[0] == 0
    , wave[a] == 0
    } /.dsol //ComplexToTrig  //Simplify
:[font = input; noPageBreak; preserveAspect]
Reduce[eq3,{C[1],C[2]}]
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
Requiring  Sin[k a] = 0,  we set k=n Pi/a where n is a positive integer, and we have  C[1] == -C[2].
:[font = input; noPageBreak; preserveAspect]
wave1[x_]= ( wave[x] 
            /.dsol 
            //.{ C[2] -> -C[1], C[1] -> I c0/2, k -> n Pi/a } 
            //ComplexToTrig  
            //Simplify
          )
;[s]
3:0,0;58,1;59,0;170,-1;
2:2,12,10,Courier,1,12,0,0,0;1,15,12,Times,0,14,0,0,0;
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
We took C[1] = I c0/2 for simplicity; we'll normalize the wave functions next. 
The constant c0 follows from normalizing the wave function.  
We get an expression for the normalization constant from   
:[font = input; noPageBreak; preserveAspect]
norm = 1/Sqrt[ Integrate[wave1[x]^2 ,{x,0,a}] ] /.{Sin[2 n Pi]->0} //PowerExpand
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
The normalized wave function is    
:[font = input; noPageBreak; preserveAspect]
psi[x_,n_,a_:1] = norm wave1[x]
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; cellOutline]
c) Determine the energy eigenvalues. 
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
We verify that these are orthonormal: 
:[font = input; noPageBreak; preserveAspect]
Table[ 
   Integrate[ psi[x,n1,a] psi[x,n2,a] ,{x,0,a}] 
   ,{n1,1,4},{n2,1,4}] //TableForm
:[font = input; noPageBreak; preserveAspect]
Plot[ Table[ psi[x,n] ,{n,1,4}] //Evaluate 
   ,{x,0,1}
   ,PlotStyle->{{RGBColor[1,0,0]}
               ,{RGBColor[0,1,0]}
               ,{RGBColor[0,0,1]}
               ,{RGBColor[0,1,1]}
               }
   ];
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
           The energy levels follow from  the expression for k, or from applying the hamiltonian operator to the wave function  and setting it equal to energy times the wave function:
:[font = input; noPageBreak; preserveAspect]
eq4= energy  psi[x,n,a] == -hbar^2/(2 m) D[ psi[x,n,a],{x,2}]
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect]
 The energy eigenvalues follow by solving for the energy 
:[font = input; noPageBreak; preserveAspect]
eSol = Solve[eq4,energy][[1]]
:[font = text; inactive; noPageBreak; keepOnOnePage; preserveAspect; endGroup; endGroup; endGroup]
 The energy levels differ by a factor of n^2.
:[font = section; inactive; noPageBreak; preserveAspect; cellOutline; backColorRed = 58981; backColorGreen = 58981; backColorBlue = 58981]
Optional Problem 3: Particle bound in an Finite Potential Well 
:[font = section; inactive; noPageBreak; preserveAspect; cellOutline; backColorRed = 58981; backColorGreen = 58981; backColorBlue = 58981]
Optional Problem 4: Modify #1 for the case of a single step potential
^*)