%%%
tic
M=1.0; %%mass
w=1;    %%trap
g=10; %%interaction
dt=0.04; 
t=0:dt:20000;
bound=10;
dx=bound/200;
x=-bound:dx:bound;

nu=0;
accuracy=1.0e-7;
Nx=length(x);

%%%initial wave function
phi0=rand(1,Nx)+1i*rand(1,Nx)-0.5-0.5*1i;
%phi0=exp(rand(1,Nx)+1i*rand(1,Nx)-0.5-0.5*1i);
nor=sum(abs(phi0).^2)*dx;
phi0=phi0/sqrt(nor); 
phi_initial=phi0;

dkx=2*pi/((Nx-1)*dx);
kx=-pi/dx:dkx:pi/dx;%Define k-space grid


chemical_mu=zeros(length(t),1);
%%
Ek=(kx).^2/(2*M);
Uk=exp(-0.5.*Ek.*dt);
Vtrap=0.5*M*x.^2*w.^2;
%%
for it=0:length(t)
    if mod(it,1000)==0
        it
        toc
    end
    Vint=g*abs(phi0).^2;
    V=Vint+Vtrap;
    U=exp(-0.5*V*dt);
   
    phi1=U.*phi0;      %%%第一步演化
    phi1k=fftshift(fft(phi1)); %%%转到k空间
    phi2k=Uk.*phi1k;      %%%第二步演化
    phi2=ifft(ifftshift(phi2k));    %%%转到实空间
    phi3=U.*phi2;      %%%第三步演化
    
    mu=-log(abs(sum(conj(phi0).*phi3)*dx))/dt; %%%化学势
    nor=dx*sum(abs(phi3).^2);
    phi3=phi3/(sqrt(nor));
   
    if abs(mu-nu)<accuracy   %%%%判断精度
        'mu_nu'
        break
    end
    nu=mu;
    chemical_mu(it+1)=mu;
    phi0=phi3;   
end	
%%

fig=figure;
set(gcf,'Position',[582 455 562 442]);
plot(x,abs(phi3).^2,'linewidth',1.6);
xl=xlabel( ['$' 'x' '$']);
yl=ylabel( ['$' '|\psi(x)|^2' '$']);
%tl=title(['$' 't=' '$',num2str(t(it+1))]);
%    tl=title(['$','t=','$',num2str(t(it))]);
set([xl,yl],'interpreter','latex');
set(gca,'fontsize',20,'linewidth',1.2)

figname=['gs_g_',num2str(g),'.fig'];
saveas(fig,figname);
%close(fig)

%%
