% thr8.m   May 1998


% This program is used in Section 6 of 
% 'Computation of Thresholds for Schrodinger Operators'
% The program enables one to find values of s for which one
% has zero energy thresholds for
% a self-adjoint operator with a coupling constant s.

%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,N);
  A2=sparse(2:N,1:N-1,-h^2*ones(1,N-1),N,N);
  Neum=sparse([N],[N],[h^2],N,N);
  % for Dirichlet boundary conditions at N use the following
 B1=A1+A2+A2';
  % for Neumann boundary conditions at N use the following
  % B1=A1+A2+A2'-Neum;
  

 
   %define the conditioning matrix T
  w=zeros(1,N);
   for r=1: N
	  x=r/h;
	 w(1,r)=1+x;
  end;
   T=sparse(1:N,1:N,w,N,N);
  
 	  
 %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,N); 
  
  
  % State the number of values of iterations to be carried out
  % in the bisection method
 M=40;
  %define the output sequence giving the smallest eigenvalue
  output=zeros(M,1);

  
  
  
  
  
  
  
  
% we start the computation which is repeated for a sequence of values of s
% between the following two values
smin=4;
smax=5;
 
 for m=1:M
	 %set the  next value of s to be tested in the bisection method
	 s=(smin+smax)/2;  
	 %set up the s-dependent matrix
	B=B1+s*B2;
%  eigg=first eigenvalue of A or of a preconditioned version of A;
% specify the matrix, then the number of eigenvalues required, then
% the value about which the chosen eigenvalues are to cluster.
	 eigg=eigs(T*B*T,1,-0.1)
  % set the values of smin and smax for the next iteration
  if eigg>0
	  smin=smin+(smax-smin)/2;
  else
	  smax=smax-(smax-smin)/2;
  end;
   %the following is to enable me to follow the progress of the computation
    disp([s])
	output(m)=eigg;
end;


s
output