clear all
close all
clc
%数值求解含时薛定谔方程，研究波函数在遇到势垒后的行为

NPTS=220;           %格点数
sigma=10;k0=0.6;    %Gauss波包参数
x0=40;              %Gauss波包参数
time=400;           %演化时间,可以用更长演化时间观测波包运动
N=NPTS-1;            %空间被离散成NPTS-1段
v(1,NPTS)=0;        %格点上的势函数V
v(1,105:115)=-12.0;   %势垒，可改变势垒高度，观察波函数穿越势垒的情况
%注意在下面的矩阵中，空间步长h=1，时间步长delta_t=1,
%系数矩阵求法见ppt中标注为(2)的公式
%A=diag(-2+2i-v(2:N),0)+diag(ones(N-2,1),1)+diag(ones(N-2,1),-1);  %系数矩阵
%使用diag函数可以使程序大大简化
A=zeros(N-1,N-1);
A(1,1)=-2+2i-v(1,2);
A(1,2)=1;
A(N-1,N-1)=-2+2i-v(1,N);
A(N-1,N-2)=1;
for n=2:N-2
    A(n,n)=-2+2i-v(1,n+1);
    A(n,n+1)=1;
    A(n,n-1)=1;
end
PHI0(1)=0;PHI0(NPTS)=0;  %边界条件
for x=2:NPTS-1;
PHI0(x)=exp(k0*i*x)*exp(-(x-x0)^2*log10(2)/sigma^2); %初始波函数，Gauss波包
end
PHI=zeros(NPTS,time);
PHI(:,1)=PHI0';
for y=2:time
CHI(2:NPTS-1,y)=4i*(A\PHI(2:NPTS-1,y-1));  %见公式(2),步长h和时间间隔delta_t都为1
PHI(2:NPTS-1,y)=CHI(2:NPTS-1,y)-PHI(2:NPTS-1,y-1);%见公式(1)
end
mo=moviein(time);
for j=1:time
plot([1:NPTS],abs(PHI(:,j)).^2/norm(PHI(:,j)).^2,'r',[1:NPTS],v,'b'); 
%画波函数时我们做了归一化处理，norm为取模运算
axis([0,220,0,0.1]);
mo(j)=getframe; %对当前图形拍照后产生的数据向量依次存放于画面构架数组中
end
movie(mo,1,50)  %以每秒50帧的速度，重复播放1次

figure;plot([1:NPTS],abs(PHI(:,1)).^2/norm(PHI(:,1)).^2,'r',[1:NPTS],v,'b')
%初始波包