usecanvas(%T);
clear
dx=0.1;
x=[-20:dx:20];
eye=sqrt(-1);
ksi=2.0;
psi=exp(-x.*x/2-eye*x*ksi)';

[n,m]=size(x);

dt=0.05;

v=x.*x*0.0;
for i=1:m
	if((x(i)>4.0)&(x(i)<5.0)) 
	  v(i)=-15.0;
	end
end

for i=2:m-1
h1(i-1,i-1)=1-dt*(v(i)+1/dx/dx)*eye;
h1(i,i+1)=dt/dx/dx*eye/2.0;
h1(i,i-1)=dt/dx/dx*eye/2.0;
end
h1(m,m)=h1(m-1,m-1);
h1(m,m-1)=h1(m-1,m-2);
for i=2:m-1
h2p(i-1,i-1)=1+dt*(v(i)+1/dx/dx)*eye;
h2p(i,i+1)=-dt/dx/dx*eye/2.0;
h2p(i,i-1)=-dt/dx/dx*eye/2.0;
end
h2p(m,m)=h2p(m-1,m-1);
h2p(m,m-1)=h2p(m-1,m-2);
h2=inv(h2p);
h1b=h2*h1;
// f=gcf();
// f.pixmap="on"
// a=get("current_axes");
// a.data_bounds=[-20,0,0,1,1,1];
i=0;
psip=psi;
for t=0:dt:5
   psip=h1b*psip;
   i=i+1;
   for j=1:m
    psit(i,j)=abs(psip(j)).*abs(psip(j));
  end  
    tt(i)=t;
end
plot3d(x',tt',psit');




