%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Integrating GPE by the method of Steepest Descent.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;

dt =.01;
g  =100;														%self-energy-ratio
R  =[-1:.05:1]*1.5;					Z  =R;
[r,z]=meshgrid(R,Z);
Psi=exp(-r.^2-z.^2);		    	Psi=Psi/norm(Psi);			%trial-Wavefunction
V  =r.^2+sqrt(8)*z.^2;											%potential
%V  =2.5*(z.^4-2*z.^2+1) + 5*r.^2;

clf;
figure(1),surf(R,Z,V)
title('\fontsize{18}External Potential');
set(gca,'FontSize',16,'LineWidth',1.2);
xlabel('\fontsize{18}{\itr/a_{ho}}');
ylabel('\fontsize{18}{\itz/a_{ho}}');
zlabel('\fontsize{18}{\itV/h\omega}');

for k=1:200
	[ gradPsi,Dumm]=gradient(Psi);
	[nablaPsi,Dumm]=gradient(gradPsi);
	H=-nablaPsi + V + g*Psi.^2;
	Psi=Psi.*(1-H*dt);
	Psi=Psi/norm(Psi);

	figure(2),surf(R,Z,Psi)
	title('\fontsize{18}BEC-Wavefunction');
	set(gca,'zLim',[0 .1],'FontSize',16,'LineWidth',1.2);
	xlabel('\fontsize{18}{\itr/a_{ho}}');
	ylabel('\fontsize{18}{\itz/a_{ho}}');
	zlabel('\fontsize{18}{\it\Psi/a_{ho}^{3/2}}');
	pause(.1);
end;
