Analytical
Inside the box, the wavefunction is:
\begin{equation} \frac{\hbar^2}{2m} \frac{d^2 \psi(x)}{dx^2} = E \psi(x) \iff \frac{d^2 \psi(x)}{dx^2} = k^2 \psi(x) \end{equation} where $k = \frac{\sqrt{2mE}}{\hbar}$. Outside the box: \begin{equation} -\frac{\hbar^2}{m} \frac{\partial^2 \psi(x)}{\partial x^2} = (E - V_0) \psi(x) \iff \frac{\partial^2 \psi(x)}{\partial x^2} = \alpha^2 \psi(x) \end{equation} where $\alpha = \frac{\sqrt{2m(V_0 - E)}}{\hbar}$. The solutions for the Schrödinger equation for a particle in a finite potential well are:
\begin{equation} \psi(x) = \begin{cases} Ce^{\alpha x}, \text{if } x < -\frac{L}{2}\\ A\sin(kx) + B\cos(kx), \text{if } -\frac{L}{2} \leq x \leq \frac{L}{2}\\ Fe^{-\alpha x} , \text{if } x > \frac{L}{2} \end{cases} \end{equation}
Next, $\psi(x)$ must be continuous and differentiable, which means that the value of the functions and their derivatives must match at the boundaries, $-\frac{L}{2}$ and $\frac{L}{2}$. Depending on whether the solution is symmetric, $\psi(-x) = \psi(x)$, or anti-symmetric, $\psi(-x) = -\psi(x)$, $A = 0$, $C = F$ or $B = 0$. $C = -F$, respectively. For the symmetric case:
\begin{equation} \begin{array}{|l@{}} Fe^{-\alpha L/2} = B\cos(k L/2)\\ -\alpha F e^{-\alpha L/2} = -\alpha B\cos(k L/2) \end{array} \end{equation}
taking the ratio of the above two to give us $\alpha = k \tan(k L/2)$. Similarly, from the anti-symmetric case we get $\alpha = -k \cot(k L/2)$.Both $\alpha$ and $k$ depend on the energy, $E$, making the above equations transcendental. To solve them we make the following substitutions:
\begin{equation} \begin{array}{|l@{}} u = \frac{L}{2} \alpha \\ v = \frac{L}{2} k = \frac{L}{2}\frac{\sqrt{2mE}}{\hbar}\quad (1) \end{array} \end{equation}
Additionally, $u^2 = u_{0}^{2} - v^2$ where $u_0 = \frac{m L^{2} V_0}{2 \hbar^2}$. This leads to:
\begin{equation} \sqrt{u_{0}^{2} - v^2} = \begin{cases} v \tan(v)\\ -v \cot(v) \end{cases} \end{equation}
To calculate the above I use the following Octave code:
clc;
clear;
close all;
x = 0 :0.1:5; % coordinates interval.
L = 1; % well width.
U_0 = 40; % well depth.
hbar = 1; % Planck's constant.
m = 1; % particle mass.
u_0sq = (m*(L^2)*U_0 ) / (2 * hbar^2);
lhs = @(x) sqrt(u_0sq - x.^2);
rhs_s = @(x) x.*tan(x);
rhs_a = @(x) -x.*cot(x);
plot(x, lhs(x), 'r', 'LineWidth', 2, ...,
x, rhs_s(x), 'b', 'LineWidth', 2, ...,
x, rhs_a(x), 'g', 'LineWidth', 2)
xlabel('v')
axis([0 5 0 5])
From $(1)$ and the three intesection values1 for $v$ from the below plot, the results for the energies are:
3.2724
12.8576
27.7294
Numerical
Now, for the numerical solution I use the following:
% Parameters for solving problem in the interval -L < x < L.
L = 5; % Interval Length.
N = 1000; % No of points.
x = linspace(-L, L, N)'; % Coordinate vector
dx = x(2) - x(1); % Coordinate step
% Finite square well of width 2w and depth given by U_0.
U_0 = 40; w = L / 10;
U = -U_0*(heaviside(x + w) - heaviside(x - w));
% Discretized (three-point finite-difference) representation of the Laplacian.
e = ones(N, 1);
Lap = spdiags([e -2*e e], [-1 0 1], N, N)
% Total Hamiltonian
hbar = 1;
m = 1;
H = -1/2*(hbar^2 / m)*Lap + spdiags(U, 0, N, N); % T = H + U.
% Find (first 3) (smalest algebraic) eigenvalues.
nmodes = 3; options.disp = 0;
[V, E] = eigs(H, nmodes, 'sa' , options);
[E, ind] = sort(diag(E));
E
The results for the energies are:
36.72
27.14
12.27
I'm expecting some error but in this case, the energies differ in an order of magnitude.
Clearly, I am doing something wrong but I don't know whether the error is in the analytical or numerical solution. What am I doing wrong here? Any hint will be appreciated.
1. The vertical lines are discontinuity points.
