clear;
% TDSE, doublewell, normal mode/energy calculation, linear mode
% calculations -- animation of tunneling
% Yael Elmatad, March 10, 2004
%(Generates Sparse nxn matrix of the form : 
%[ 2 -1  0]
%[-1  2 -1]
%[ 0 -1  2]
nn = 201;
[n,x1,bh,k,u,E,V,D] = efuncs (nn);
Vx = [u.*x1.^4 - k.*x1.^2 + bh];
ne = 10;
y = V(1:n, 1:ne);

% animation... reuse above code with t increasing in a loop... 
% lc = linear combination
% b1, b3 = coefficients (normalized)
figure;
for t=0:0.1:100
  ss1 = 5; %normal mode number - stationary state
  ss2 = 6; %second normal mode number - stationary state 2
  b1 = (1./2).^(1./2); %coefficient of first state
  b2 = (1 - b1.^2).^(1./2); %coefficient of second state 
  Elc = (b1.^2)*(D(ss1,ss1)) + (b2.^2)*(D(ss2,ss2));  %energy of new state
  lc = b1.*V(1:n, ss1).*exp(-i.*D(ss1,ss1).*t)+ b2.*V(1:n,ss2).*exp(-i.*D(ss2,ss2).*t); %linear combination
  plot (x1, Vx); %plots doublewell potential
  hold on; plot(x1, 200*[real(lc) imag(lc)] + Elc); hold off;
  axis([-1 1 0 750]);   % need to fix range otherwise will jump around
  title('Time dependant double well potential');
  xlabel('x'); ylabel('u(x,t)'); text(1,1, sprintf('t=%g', t));
  drawnow;
  pause(.1);
end