Since the problem is rotationally symmetric rather than square-symmetric, it makes more sense to use a coordinate system that reflects the symmetry in the problem. Hence, I'd like to solve this problem is spherical coordinates (r'--).
In spherical coordinates Schrödinger's equation reads:
Inside the well, where the potential is zero, U0 would be set to zero.
The stuff in the square brackets (with the overall minus sign) is the (dimensionless) angular momentum squared operator: L'2. The eigenfunctions of this operator are degenerate, so there is a choice to make in defining the basis states. In physics folks conventionally use the Ylm, which are simultaneously eigenfunctions of both L'2 and L'z:
L'2 Ylm = l(l+1) Ylm
L'z Ylm = m Ylm
I direct your attention to the H-atom and SHM, where we have discussed a bit about these Ylm.
Taking a as our length scale and 2/(2ma2) as our energy scale, we gain the dimensionless form
Seeking a solution where the dependence on r' factors from the dependence on angles ("separation of variables", =R(r') Ylm(,)), we find:
We will find it convenient to divide the above equation by E'-U0', but we will need to pay attention as to the sign of (E'-U0'). When (E'-U0')>0 we are in a classically allowed region, and we expect the wavefunction to oscillate with a wavelength that depends on the particle's momentum... sin(kr) like solutions are to be expected. When (E'-U0')<0 we are in a classically disallowed region, and we expect the wavefunction to exponentially decay...exp(-r) like solutions are to be expected. We'll handle both cases at the same time, sometimes by using ± signs where the top case (+ here) refers to a classically allowed region and the bottom case (- here) refers to a classically disallowed region. We substitute:
in our r' differential equation, producing a differential equation which I'll write in three equivalent forms:
As usual we'll have to work to solve the radial differential equation! Start by seeing how the the equation must work for large . If in the last equation we ignore the (very small) -2 term, we get the familiar oscillator differential equation. For E>U0 we find:
That is we see the expected oscillations as a function of r', but for large r' with ever smaller amplitude, because of the 1/r factor.
For E<U0 we find:
That is we see the exponential behavior as a function of r' expected in a classically disallowed region.
For small , the -2 term dominates and we can ignore the right hand side in comparison. If we try a power-law solution n
we find n=l.
Factoring out all the required behavior for small , we hope to find a simple power series for R:
Plugging this in to the first R equation:
The result is a two term recursion relation (i.e., the result has just two as so, for example, given a0 we can calculate a2, from which we can calculate a4, etc. until we're done.) Note that the the recursion relation connects even k to even k and since a0 may not be zero (as the leading term l was needed to make the small differential equation work). Thus the series will have only even terms and we write our even k=2i, where i is an integer.
The series does not terminate; the result (for E>U0) can be recognized as the spherical Bessel function jl from the hypergeometric formula (0F1).
As shown above, spherical Bessel functions are related to normal Bessel functions with half-integral index.
As we expect that a second order linear differential equation should have two independent solutions, we're only half way there. jl is regular at =0, the other spherical Bessel function: yl (also known as nl) explodes at =0, and hence is not part of normalizable wavefunctions. Here are the basic large and small argument behaviors of jl and yl.
As we showed above, for large we expect oscillations like exp[±i]/. For small we expect l or -l-1, exactly as shown above. Here are some plots:
y and j are real functions and hence oscillate like sin and cos. If we want, we can put y and j together to produce something that for large oscillates like exp[±i]/. The results are related to Hankel functions and are also called spherical Bessel function of the third kind:
For E<U0 we've found the exponentially rising modified spherical Bessel function related to Il+½:
As above there is a second solution that is exponentially falling related to Kl+½:
The above expression may make spherical Bessel functions look as complex as the generic Bessel functions, but in fact, the spherical Bessel functions can be expressed in terms of well-known functions (whereas a generic Bessel function can not be so easily expressed):
jl(k'r')=0 for r'=1
That is k' must be a zero of jl. The energy of the wavefunction can then be calculated from
Here is the overall solution:
Because the angular dependence is "simple" we can usefully plot the wavefunction just as a function of r'. Here are "stacked wavefunction" plots for l=0,1,2:
The red line is the classical "effective potential": l(l+1)/r'2 (which is the "centrifugal potential") plus the infinite square well. Note that the centrifugal barrier excludes non-zero l wavefunctions from the origin.
You should check each of the above roots on the above plots of jj.
Here is what the set of energy levels looks like: