# Discrepancies between numerical and analytical solution for particle in a finite potential well?

## 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.

I think that your potential for the numerical case is wrong. The potential should be a big positive number, so the solution tends to zero outside the well when the value of the potential increases.

If you check the solutions in Wikipedia for the first three (bounded) states they should be

$$E_n = \frac{2\hbar^2 v_n^2}{m L^2}\, ,$$

with $v_1 = 1.28$, $v_2 = 2.54$ and $v_3=3.73$, or just

\begin{align} &E_1 = 3.2768\\ &E_2 = 12.9032\\ &E_3 = 27.8258 \end{align}

Following is your code with corrections, notice that you also miss the step size of the finite difference.

% 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*(1 - 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)/dx^2;

% Total Hamiltonian
hbar = 1;
m = 1;
H = -0.5*(hbar^2 / m)*Lap + spdiags(U, 0, N, N);  % T = H + U

% Find (first 3) (smalest algebraic) eigenvalues.
nmodes = 6;
options.disp = 0;
[V, E] = eigs(H, nmodes, 'sm' , options);
[E, ind] = sort(diag(E));
disp(E)


The following is the answer given by the code

 3.2724
12.8576
27.7294
40.1985
40.2454
40.7874