%****************************************************************
% Program 4:本程序完全从数值的角度出发,把求解薛定谔方程的过程看成是一个
% 求解哈密顿数值矩阵的特征值以及对应特征向量的一个过程
%****************************************************************
% 确定参数|求解薛定谔方程中使用到的一些参数:
hbar=1;m=1; %在原子单位体系里(hbar= 1,质量m= 1)
L = 5; % 区间长度为2L:-L < x < L
N = 1000;x = linspace(-L,L,N)'; % 总共在区间内均匀地取N个点
dx = x(2) - x(1); % 步长大小
% 选定U-x关系|在这里可以任意以多项式的形式进行组合 (这里首先假定的是一维谐振子的情况)
U = 1/2*100*x.^2; %一维谐振子满足的U-x关系
%U = 1/2*100*x.^2+1/2*x.^(4); %经过多项式组合的U-x关系(可以把上面U-x注释掉，替换成这个U-x关系)
%读者也可以自己以多项式的形式组合
%写出哈密顿的数值表达式
e = ones(N,1); Lap = spdiags([e -2*e e],[-1 0 1],N,N)/dx^2; %二阶导数的数值表达形式
H = -1/2*(hbar^2/m)*Lap + spdiags(U.*e,0,N,N); % 完整哈密顿矩阵
%计算能量本征值与能量本征函数
nmodes = 3; options.disp = 0;
[V,E] = eigs(H,nmodes,'sa',options); %利用求解特征值与特征向量的方式求出本征值与本征函数
[E,ind] = sort(diag(E));% 将E按数值大小从低到高排序
V = V(:,ind); %取出前三级的本证函数
Usc = U*max(abs(V(:)))/max(abs(U)); %为了便于观察，将U的范围按比例缩小U-x关系不变
%展示结果(如果读者认为展现的图片看不清其中的曲线，请调整显示的曲线数目)
plot(x,V,x,Usc,'--k'); %展现前三级在一定的U-x关系下本征函数V-x的关系
lgnd_str = [repmat('E=',nmodes,1),num2str(E)]; %添加标注
legend(lgnd_str)
shg %显示图形窗口
%%读者可以取消下面几行程序的注释,验证求得的能量本征函数是否满足归一化
% SumV = 0;
% for numV=1:1000
% SumV = SumV+V(numV,1)^2;
% end
% SumV
