
function psiright=double_well(E)
% this is a double well potential,  extending from x=-a-b to x=a+b.  The
% potential is zero, except in the region -a to +a, where it is Vo

global flags a b mass Vo

% if flags is +1, then symmetric state
% if flags is -1, then antisymmetric state
  
 % wavevectors
 if E>Vo
    % classically allowed region
      kiii=sqrt(2*mass*(E-Vo));
      if flags==1 % symmetric state, wavefunction in region III is cos(kiii x)
         bvec=[cos(kiii*b/2) -kiii*sin(kiii*b/2)]';
      elseif flags==-1 % antisymmetric state, wavefunction in region III is sin(kiii x)
         bvec=[sin(kiii*b/2) kiii*cos(kiii*b/2)]';
      end
 else
    % classically forbidden region     
      kappa=sqrt(2*mass*(Vo-E));
      if flags==1 % symmetric state, wavefunction in region III is cosh(kiii x)
         bvec=[cosh(kappa*b/2) kappa*sinn(kappa*b/2)]';
      elseif flags==-1 % antisymmetric state, wavefunction in region III is sin(kappa x)
         bvec=[sinh(kappa*b/2) kappa*cosh(kappa*b/2)]';
      end
end
 k2=sqrt(2*mass*E);
     
 Cmat = [cos(k2*b/2) sin(k2*b/2); k2*[-sin(k2*b/2) cos(k2*b/2)]];

 % determine coefficients, A and B

 avec=Cmat\bvec;
 
 % wavefunction at the right-hand boundary

 psiright=avec(1).*cos(k2*(a+b/2))+avec(2).*sin(k2*(a+b/2));
 
 
 