Numerical solutions to 1D Schrodinger equation suggest degenerate energy eigenvalues, though this is supposedly disallowed

Question

I am working on solving the time-independent Schrodinger equation using the method of finite differences. This approach has been discussed previously on this site (here, for instance). My code is written in Mathematica and is able to reproduce the energy eigenvalues and wavefunctions for a 1D particle-in-a-box and the harmonic oscillator.

When I take the particle-in-a-box potential and add in a large barrier in the middle third of the box, such that $$ V(x) = \begin{cases} 0, & 0 < x < \frac{1}{3} \text{ and } \frac{2}{3} < x < 1 \\ 2000, & 1/3 \leq x \leq 2/3 \\ \infty, & \text{otherwise} \end{cases}$$

I get an interesting result. The solutions appear to be particle-in-a-box wavefunctions that are constrained to exist in either the left or right partitions (see picture below). Now, I know that quantum tunneling causes some probability density to leak through, but I don't believe that this is important. The oddity is that the wavefunctions come in pairs, and their energy eigenvalues appear to be pair-wise degenerate.

According to problem 2.44 in the third edition of Quantum Mechanics (Griffiths and Schroeter), the 1-D Schrodinger equation does not permit degenerate eigenvalues for linearly-independent solutions. There are a few assumptions that go into the derivation, one of which is that the solutions approach 0 as $x \rightarrow \infty$. This is certainly the case here, and $\psi_1$ and $\psi_2$ are definitely linear-independent.

With all of that said, where is the issue? Is the theorem wrong? Am I not accounting for an assumption that is violated in this problem? Or does it have to do with some voodoo relating to the indistinguishability of electrons or entanglement? Or perhaps the solutions are linearly-dependent in phase space, and Griffiths over-simplified the wording of the theorem?

Either way, I'd love some insight.

EDIT:

I may have found the reason. When I force Mathematica to compute the eigenvalues using infinite precision, then I get symmetric or antisymmetric wavefunctions, as expected. I only tried this for small matrix sizes (50), and will try to confirm with a larger size.

Here's my Mathematica code. I apologize if it is messy. I am a python programmer and did just enough to get this thing to run:

(* Matrix size *)
n = 100;

(* Range = 0 < x < a *)
a = 1;

(* Calculate \[CapitalDelta]x *)
\[CapitalDelta]x = a / (n + 1);

(* Specify \[Nu](x) = V / \[Lambda] where V = potential function and \
\[Lambda] = \[HBar]^2/(2m(\[CapitalDelta]x)^2) *)
\[Nu][x_] = If[x <= a/3 || x >= 2 a/3, 0, 2000];

(* Create and populate array. *)
data = Array[0 &, {n, n}];
For[i = 0, i < n, i++,
    For[j = 0, j < n, j++,
        If [i == j, 
   data[[i + 1, j + 1]] = 2 + \[Nu][(j + 1) * \[CapitalDelta]x ]  , 0];
        If [j == i - 1 || j == i + 1, data[[i + 1, j + 1]] = -1, 0];
    ]
 ]
N[data] // MatrixForm;

(* Find all eigenvectors and eigenvalues. Sort them by ascending \
energy values *)
eigensystem = Eigensystem[data, -n];
eigensystem = Transpose@SortBy[Transpose[eigensystem], First];


(* Define state. k = 1 represents ground state *)
k = 3;

(* Extract energy value. *)
energy = \[HBar]^2/(2 m (\[CapitalDelta]x)^2)*eigensystem[[1]][[k]]
N[energy]

(* Extract eigenvector of choice (and append 0 in beginning and end *)


eigenvector = Insert[eigensystem[[2]][[k]], 0.0, {{1}, {n + 1}}];

(* Plot potential, wavefunction, and probability density *)
plot1 = Thread[{Subdivide[a, n + 1], eigenvector}];
plot2 = Thread[{Subdivide[a, n + 1], eigenvector^2}];
GraphicsRow[{Plot[\[Nu][x], {x, 0, a}, PlotRange -> All], 
  ListLinePlot[plot1, PlotRange -> All], 
  ListLinePlot[plot2, PlotRange -> All]}, ImageSize -> 1000]


(* Plot a range of wavefunctions and probability densities *)
kstart = 1;
kend = 4;
GraphicsRow[
 Table[ListLinePlot[
   Thread[{Subdivide[a, n + 1], 
     Insert[eigensystem[[2]][[i]], 0.0, {{1}, {n + 1}}]}], 
   PlotRange -> All], {i, kstart, kend}], ImageSize -> 1500]
GraphicsRow[
 Table[ListLinePlot[
   Thread[{Subdivide[a, n + 1], 
     Insert[eigensystem[[2]][[i]]^2, 0.0, {{1}, {n + 1}}]}], 
   PlotRange -> All], {i, kstart, kend}], ImageSize -> 1500]
```

Answer 1

Running your Mathematica code with using Eigensystem[N[data], -k] results in your original results while with precise values the eigenstates are symmetric/antisymmetric (but any large size will take a much longer time for Mathematica to compute).

This is all an issue in what Mathematica rounds as an eigenvector: you can mathematically show the approximation matrix $A$ (data) used is nondegenerate and does not have any eigenvectors that aren't symmetric or antisymmetric by its inherent structure. If one assumes $\psi$ is an eigenvector of $A$ with eigenvalue $\lambda$, then given a fixed value of $\psi_0$ the necessary value of $\psi_1$ can be uniquely determined ($\lambda \psi_0 = 2 \psi_0 - \psi_1$ turns into $\psi_1 = (2 - \lambda)\psi_0$), and then $\psi_2$ can be uniquely found given $\psi_0$ and $\psi_1$ through re-arranging $\lambda \psi_1 = -\psi_0 + 2\psi_1 - \psi_2$ from there and so on since the effect of $A$ on an element of $\psi$ always relies on its nearest neighbors, even within the middle potential well. This implies vectors corresponding to the same eigenvalue of $A$ can only differ multiplicatively, and thus the eigenstates are nondegenerate. Since running this logic backwards from the last element works exactly the same (since the potential is centered), and any multiplication on a starting element exactly carries over to all generated elements, $\psi_0$ and the last value $\psi_f$ can only differ in sign or they would not result in each other for both inductions. Thus, for this discrete approximate Hamiltonian the energy levels are not degenerate and are either symmetric or antisymmetric over the center.