%  thr13.m    May 1998


% This program is used in Section 6 of
% Computation of Thresholds for Schrodinger Operators".
% The program enables one to find the eigenvalues of the modified operator
% THT when the coupling constant is set at the value for which there is a 
% zero energy threshold. 

% The matrix includes the tail correction term.


%semi-colons stop output into the command window
close('all')
%clear out all the variables
clear variables;



%  Define the basic constants used.
% h=number of divisions per unit length
h=20
%b=length of interval
b=10
% N=dimension of space
N = h*b



%define the kinetic energy matrix
u=ones(N,1);
v=ones(N-1,1);
% SPARSE MATRICES
%the first two vectors below give the positions of the 
% non-zero entries and the third vector gives their values.
% the last two numbers give the size of the matrix
A1=sparse(1:N,1:N,2*(h^2)*ones(1,N),N+1,N+1);
A2=sparse(2:N,1:N-1,-(h^2)*ones(1,N-1),N+1,N+1);
B1=A1+A2+A2';



%define the conditioning matrix T
w=zeros(1,N+1);
for r=1: N+1
	x=r/h;
	w(1,r)=1+x;
end;
T=sparse(1:(N+1),1:(N+1),w,N+1,N+1);



%define the potential sequence
pot=zeros(1,N);
for r=1: N
	x=r/h;
	pot(1,r)=cos(pi*x)*exp(-x^2/2);
end;
%define the potential energy matrix
B2=sparse(1:N,1:N,pot,N+1,N+1); 








% determine the value of the normalisation constant for the tail vector
cc=0;
for r=1:(N+h)
	cc=cc+r^(-2);
end;
ccc=(pi^2)/6-cc;
c=ccc^(-1/2)/h;



% find the terms of the matrix associated with the tail vector
C2=sparse([N+1],[N+1],[c^2*h^2],N+1,N+1);
entry=-c*h^2*w(1,N)
C3=sparse([N],[N+1],[entry],N+1,N+1);
C4=sparse([N+1],[N],[entry],N+1,N+1);



% set up the s-dependent matrix C
% State the initial values of s and then
% construct an iterative loop which determines the 
% threshold value of s more and more accurately
smin=4.25
smax=4.26
for t=1:40
	s=(smin+smax)/2
B=B1+s*B2;
C1=T*B*T;
C=C1+C2+C3+C4;



%  eig12=first two eigenvalues of C;
% specify the matrix, then the number of eigenvalues required, then
% the value about which the chosen eigenvalues are to cluster.
eig12=eigs(C,2,-0.00001);
one=eig12(1)
two=eig12(2)



% we redefine smin or smax according to the sign
% of one
if one>0
	smin=s;
else
	smax=s;
end;
% we now end the loop determining the threshold s
end;



%  eigN=largest eigenvalue of C;
f=zeros(N+1,1);
f(N)=1;
for r=1:10 
 	g=C^32*f;
	f=g/norm(g);
end;
eigN=norm(C*f)/norm(f)
%finally compute the relative condition number
condition=eigN/two