The wavefunction ψ(ξ) of a particle in a box satisfies the Schrödinger equation

-ψ''=εψ ,

where the prime denotes the derivative over the dimensionless variable ξ=x/L, x is the coordinate of the particle and L is the size of the box, ε is the dimensionless energy. The actual energy is E=(ℏ²ε)/(2mL²).

The boundary condition is ψ(0)=ψ(1)=0.

Numerically one can represent the wave-function on a grid

ξ_i = ih, h=1/N, i=0...N,

where the wave-function is represented by a vector

ψ_i = ψ(x_i), i=0...N.

The second derivative operation can then be approximated using the three point finite difference formula,

ψ''_i = (ψ_i-1-2ψ_i+ψ_i+1)/h² .

Now the Schrödinger equation becomes a matrix eigenvalue equation

∑_j=1,N-1 H_ij ψ_j = ε ψ_i, i=1...N-1,

where—since the wave-function at the end-points is zero—the summation goes over 1...N-1 and vanishing components ψ₀ and ψ_N have been omitted.

The matrix H in the three-point approximation is given as (taking into account the zero boundary condition and omitting the first and the last component of the wave-function vector)

	(	-2	1	0	0	0	)
	(	1	-2	1	0	0	)
	(	0	1	-2	1	0	)
-h²H =	(	...	...	...	...	...	)
	(	0	0	1	-2	1	)
	(	0	0	0	1	-2	)

Solving the matrix eigenvalue equation amounts to matrix diagonalization.

1. Build the Hamiltonian matrix h²H as a gsl_matrix structure. Hint:

   int n=20;
   gsl_matrix *H = gsl_matrix_calloc(n,n);
   for(int i=0;i<n-1;i++){
     gsl_matrix_set(H,i,i,-2);
     gsl_matrix_set(H,i,i+1,1);
     gsl_matrix_set(H,i+1,i,1);
     }
   gsl_matrix_set(H,n-1,n-1,-2);
   gsl_matrix_scale(H,-1);
   
   

2. Diagonalize the matrix using GSL subroutines for diagonalization of real symmetric matrices. Hint:

   gsl_eigen_symmv_workspace *w =  gsl_eigen_symmv_alloc (n);
   gsl_vector *eval = gsl_vector_alloc(n);
   gsl_matrix *evec = gsl_matrix_calloc(n,n);
   gsl_eigen_symmv(H,eval,evec,w);
   

3. Sort the eigenvalues and eigenvectors.

   gsl_eigen_symmv_sort(eval,evec,GSL_EIGEN_SORT_VAL_ASC);
   

4. Check that your energies are correct. Hint: the scale should be close to unity (why?):

     for(int k=0;k<n/3;k++){
       double scale=gsl_vector_get(eval,k)*(n+1)*(n+1)/pi/pi/(k+1)/(k+1);
       fprintf(stderr,"%g\n",scale);
       }
   

5. Plot the several lowest eigenfunction (and compare with the analytical result).

   for(int k=0;k<3;k++){
     printf("%i %f\n",0,0.0);
     for(int i=0;i<n;i++) printf("%i %f\n",i+1,gsl_matrix_get(evec,i,k));
     printf("%i %f\n",n+1,0.0);
     printf("\n\n");
     }
   

6. You may try to investigate how the accuracy changes with the number of grid points.
7. You may try impove the accuracy of the method by using the five-point formula for the finite difference approximation to second derivative from here.