clear;
% TDSE, tunneling time plot, fd approx
% Yael Elmatad, April 13, 2004
%(Generates Sparse nxn matrix of the form : 
%[ 2 -1  0]
%[-1  2 -1]
%[ 0 -1  2]
nn= 200; %select an even number 
[n,x1,bh,k,u,E,V,D] = efuncs (nn);
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)
  f = 1:1:50; %amount of linear combinations
  ss2 = 2.*f; %second normal mode number - stationary state 2
  ss1 = 2.*f-1; %normal mode number - stationary state
  b1 = (1./2).^(1./2); %coefficient of first state
  b2 = (1 - b1.^2).^(1./2); %coefficient os second state 
  Elc = diag((b1.^2).*(D(ss1,ss1)) + (b2.^2).*(D(ss2,ss2)));  %energy of new state
  delE = diag(D(ss2,ss2) - D(ss1,ss1)); %energy spacing
  tunt = pi./delE; %tunneling time
  Er = Elc - bh;
  figure;
  semilogy(Er, tunt); %plots doublewell potential
  title('Tunneling Time vs. <E>-bh');
  xlabel('[<E>-bh]'); ylabel('log t');
  axis ([-500 20000 .01 50000]);