function Schrodinger1DTDbnew(Flag,N,Nm);
% Schrodinger1DTDbnew('Square well',800,5)
% Schrodinger1DTDbnew('Step well',800,5)
%

% length units: A
% Time unit: 0.1 fs


E=5;   % ev
V=2*E;

L=300;
Imax=600;
dx=L/Imax;

Well_width=0.1*L;      % in meter    
x0=0.25*L;
sigma0=18*dx;

%============================================================

k0=sqrt(0.263*E);
dt=0.3;
xdt=dx*dx/dt;
azi=3.457*xdt;
azi2=2*azi;

Ic=round(0.5*Imax);
Ib1=round(0.5*Imax)
Ib2=Ib1+round(Well_width/dx)


Potential=zeros(Imax,1);
switch Flag
case 'Step well'
    Potential(Ic:end)=V;    
case 'Square well'
    Potential(Ib1:Ib2)=V;    
end


%=======================================================
% Initialize wave function

x=(1:Imax)*dx;
x=x';
Psi=exp(i*k0*(x-x0)).*exp(-((x-x0)/sigma0).^2/2);

% normalization factor (small mistake!)
%Psi2=sqrt(sum(abs(Psi).^2));
%Psi=Psi/Psi2;

% normalization factor (correct)
Psi2=sqrt(sum(abs(Psi).^2)*dx);
Psi=Psi/Psi2*dx;


kk=0:k0/1000:2*k0;
Ak=exp(-(kk-k0).^2*sigma0^2/2);
figure(1)
plot(kk,abs(Ak))

%=======================================================
T=zeros(Imax,Imax);
for I=1:Imax
    if I==1 
        T(I,Imax)=1;
        T(I,I+1)=1;
    elseif I==Imax
        T(I,I-1)=1;
        T(I,1)=1;
    else
        T(I,I-1)=1;
        T(I,I+1)=1;
    end        
    T(I,I)=-2-0.2627*dx*dx*Potential(I)+i*azi;
end
InvT=inv(T);

s=0.025;
maxPsi=max([max(abs(Psi)),s*E,s*V])*1.2;

figure(2)
clf;
line('xdata',x,'ydata',Potential*s,'color','g','linestyle','--','linewidth',2)
h1=line('xdata',x,'ydata',real(Psi),'color','r');    
h2=line('xdata',x,'ydata',imag(Psi),'color','b');    
h3=line('xdata',x,'ydata',abs(Psi).^2,'color','k');    

xlabel('x (A)'); ylabel('Potential (ev)')
box on
axis([min(x) max(x) -maxPsi maxPsi])

MM(1)=getframe;
NormC0=sum(abs(Psi).^2);

%=======================================================
%  Time marching

for n=1:N

    kai=i*InvT*azi2*Psi;
    %kai=linsolve(T,2*Psi);
    
    Psi=kai-Psi;
    
    if mod(n,Nm)==0
        NormC=sum(abs(Psi).^2)/NormC0;
        set(h1,'ydata',real(Psi))
        set(h2,'ydata',imag(Psi))
        set(h3,'ydata',abs(Psi).^2)
        title(['Time step = ', num2str(n),'     Normalization Constant=' num2str(NormC)]);
        drawnow;
        MM(end+1)=getframe;
    end
    
end

movie2avi(MM,'T6.avi')