%x^2 potential: qm result the same with gaussian packet!
%x^3 not the same anymore: make more wf narrow and they approach each other
%x^2 potential and gaussian superposition: not same anymore!
clear all
%dp/dt=-<grad V>

   xsqr_pot=2;

   n=1001;
xmin=-10;xmax=-1*xmin;
dx=(xmax-xmin)/(n-1);

cnum=20;
cposmin=0;cposmax=4;
width=0.5*4;



for ci=1:cnum
ci
	 cpos(ci)=cposmin+(cposmax-cposmin)*(ci-1)/(cnum-1);


%pot=x^2-->> grad=2*x
%pot=x^3-->> grad=3*x^2
  if (xsqr_pot==1)
classical_force(ci)=-2*cpos(ci);
else
classical_force(ci)=-3*cpos(ci)^2;
end

  for in=1:n
	   x(in)=xmin+dx*(in-1);
if (xsqr_pot==1)
	   force(in)=-2*x(in);
 else
force(in)=-3*x(in)^2;
end
%Wavefunction
wf(in)=exp(-1*(x(in)-cpos(ci))^2/(2*width^2));
%wf(in)=exp(-1*(x(in)-cpos(ci))^2/(2*width^2))+exp(-1*(x(in)+cpos(ci))^2/(2*width^2));

end
wf=wf/sqrt(sum(abs(wf).^2*dx));


%Expectation of the force at cpos
	   qm_force(ci)=0;
for in=1:n
	   qm_force(ci)=qm_force(ci)+force(in)*abs(wf(in))^2*dx;
end
qm=qm_force(ci)




figure(1);clf;
fs=20;
subplot(2,1,1);plot(x,force,'LineWidth',2);grid on;
xlabel('x','FontSize',fs);ylabel('F_{classical}','FontSize',fs);
%axis([-3 3 -30 0])
a=gca;set(a,'FontSize',fs);
subplot(2,1,2);plot(x,wf,'LineWidth',2);grid on;
xlabel('x','FontSize',fs);ylabel('\psi(x)','FontSize',fs);
a=gca;set(a,'FontSize',fs);
pause(0.2);

end  %ci


figure(2);clf;
fs=20;
subplot(1,1,1);
h1=plot(cpos,qm_force,'r','LineWidth',2);grid on;
hold on;
h2=plot(cpos,classical_force,'b--','LineWidth',2);grid on;
xlabel('position','FontSize',fs);ylabel('Force','FontSize',fs);
a=gca;set(a,'FontSize',fs);
legend([h1 h2],'Quantum','Classical');

%deviation=(classical_force-qm_force)./classical_force;
%
%subplot(2,1,2);plot(cpos,deviation,'k','LineWidth',2);grid on;
%xlabel('position','FontSize',fs);ylabel('Deviation','FontSize',fs);
%a=gca;set(a,'FontSize',fs);