Physics Simulation Forum

Dear friends,
I have written a matlab code to solve the two-dimensional schrodinger equation in a 2D Harmonic potential.
But the program gives diverging result.
Can anybody help me to find out the bug in the program. I am attaching the Matlab program.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clear all
m=200;
n=m;
T=50;
dt=0.2;
hbar=1;
me=1;
xstart=-2;
xend=2;
ystart=-2;
yend=2;
dx=(xend-xstart)/(m-1);
dy=(xend-xstart)/(m-1);
x=(xstart:dx:xend);
y=(ystart:dy:yend);
%%%%%%%Two-Dimensional_Potential%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[xx,yy]=meshgrid(x,y);
Vxy=2*(xx.^2+yy.^2);
%figure,surf(Vxy)
%shading flat
%%%%%%%%%%%%%%%%%%initial wave packet%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vecini=0.01*exp(-xx.*xx-yy.*yy);
vec=vecini;
%figure,surf(vecini)
%shading flat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dx=sparse(n-1,n-1);
Dy=sparse(n-1,n-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rx=hbar*dt/(4*me*dx*dx);
ry=hbar*dt/(4*me*dy*dy);
cv=dt/(2*hbar);
%%%%%%%%%%%%%%%%%%%Tridiagonal Matrices%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dx=diag(rx*ones(1,n))+diag((i-2*rx)*ones(1,n-1),-1)+diag(rx*ones(1,n-1),1);
Dy=diag(ry*ones(1,n))+diag((i-2*ry)*ones(1,n-1),-1)+diag(ry*ones(1,n-1),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vec(1,:)=0; %Dritchlet boundary conditions
vec(n,:)=0;
vec(:,1)=0;
vec(:,n)=0;
for tt=1:T
%starting peaceman ADI method%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%Y-direction derivative%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:length(x)-1
g(1)=0;%Dritchlet boundary conditions
g(n)=0;
for k = 2:length(y)-1
g(k)=(i+2*ry)*vec(j,k)-ry*vec(j,k+1)-ry*vec(j,k-1)+dt/2*Vxy(j,k)*vec(j,k);
end
vecn(:,j)=Dx\g';
end
vecn(:,1)=zeros(n,1);
vecn(:,n)=zeros(n,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%X-direction derivative%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 2:length(y)-1
g(1)=0;
g(n)=0;
for j = 2:length(x)-1
g(j)=(i+2*rx)*vecn(k,j)-rx*vecn(k,j+1)-rx*vecn(k,j-1)+dt/2*Vxy(k,j)*vecn(k,j);
end
vecn2(:,k-1)=Dy\g';
end
vecn2(:,1)=zeros(n,1);
vecn2(:,n)=zeros(n,1);
vec=vecn2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
surf(x,y,abs(vec.*conj(vec)))
shading flat
F(tt)=getframe()
end
movie(F)
figure,surf(x,y,Vxy),shading flat,title('2D harmonic potential')
figure,surf(x,y,vecini),shading flat,title('initial gaussian wavepacket')