Solution

The equation in question is y''(r)+V(r)y(r)=E y(r). The idea behind the matrix way of solving a differential equation is to consider the left hand side of this equation as an operator acting on the function y, i.e. L[y]=E y, and then transform the operator into matrix form with the method of finite differences. Then we have a matrix eigenvalue equation, denoted by M Y=E Y, which is easy to solve using a computer mathematics program. The eigenvectors of M, i.e. the vectors Y, are vectors with the values of the function y as their elements, i.e. (Y = (y(r_0), y(r_1), ..., y(r_n)))^T. The second derivative of y is transformed as
    y''(r_i)=(y(r_ (i + 1)) - 2y(r_i) + y(r_ (i - 1)))/h^2=(y_ (i + 1) - 2y_i + y_ (i - 1))/h^2,
so Schrödinger's equation in finite difference form looks like
    -(y_ (i + 1) - 2y_i + y_ (i - 1))/h^2 + V_iy_i = E y_i  -1/h^2y_ (i - 1) + (2/h^2 + V_i) y_i - 1/h^2y_ (i + 1) = E y_i
The idea is to form this equation at each of the points {r_i}. Then we have a set of n linear equations for the variables y_i, which can be written as M Y = E Y, and the elements of the matrix M are m_ (i, j)="the coefficient of y_j in the i:th equation". For example m_ (i, i + 1)=-1/h^2, m_ (i, i)=2/h^2+V_i. Most of the elements of M are zero. In fact, only the diagonal and the adjacent bands contain non-zero elements:
M=(

m_11 m_12 0 0 0
m_21 m_22 m_23 0 0
0 m_32 m_33 m_34 0
0 0 m_43 m_44 m_45
0 0 0 m_54 m_55
).

First, define the interval [r_a,r_b] in which to search for the solution, and the spacing of the points {r_0=r_a,r_1=r_a+h,r_2=r_a+2h,...,r_n=r_b} on this interval:

In[29]:=

h = 0.05 ;

ra = 0. ;

rb = 20. ;

n = (rb - ra)/h + 1 ;

The potential is of the form V(r)=r^2+ε r^4, with ε=0.1. Calculate its values in points {r_i}:

In[33]:=

ϵ = 0.1 ;

V[r_] := r^2 + ϵ r^4 ;

vtab = Table[V[r], {r, ra, rb, h}] ;

Then, construct the matrix M. The only non-zero elements for the ith row are m_ (i, i - 1), m_ (i, i) and m_ (i, i + 1).

In[36]:=

M = Table[Switch[i - j, -1, -1/h^2, 1, -1/h^2, 0, 2/h^2 + vtab[[i]], _, 0], {i, n}, {j, n}] ;

Now, the solutions of the equation M Y=E Y, the eigenvalues E and the corresponding eigenvectors Y can be calculated:

In[37]:=

eigvec = Eigenvectors[M] ;

eigval = Eigenvalues[M] ;

Then, we can plot the eigenfunctions, which can be obtained by means of interpolation.

In[39]:=

np = 5 ;

Plot[Evaluate[Table[Interpolation[eigvec[[n - i]]][(x - ra)/h + 1], {i, 0, np}]], {x, ra, 5}, PlotRangeAll, PlotStyleTable[Hue[i/np], {i, 0, np}], ImageSize700] ;

[Graphics:../HTMLFiles/ex10_solutions_99.gif]

Created by Mathematica  (April 10, 2007) Valid XHTML 1.1!