function [Er,params] = fnal2(u,m,params)
% Computes the eigenvalue equation of the SM
% introduced in  Fateev and Onofri, Fermilab-Pub-03/165-T
% Exploits sparse matrices, Richardson extrapolation
% and eigs(A,B,n,OPTS) to solve generalized symmetric eigenvalue problem
% Usage:
%    E = spectrum(u, m, params)
% 'params' is a structure defined as
%  params.lancz = 20;  dimension of Lanczos matrix
%  params.order = 2;   scaling order for Richardson extrapolation
%  params.side = 50;   -side/2 < x < side/2
%  params.neig = 3;    number of eigenvalues computed
%  params.maxexp = 11;   2^minexp <= dim <= 2^maxexp
%  params.minexp = 6
%  params.plot = 'y' to display eigenfunctions and potential
%
%  FERMILAB, May 2003

if(nargin<3)
  params.lancz  = 60;
  params.order  = 1;
  params.side   = 24;
  params.neig   = 10;
  params.maxexp = 14;
  params.minexp = 8;
  params.plot   = 'y';
end

lancz  = params.lancz;
order  = params.order;
side   = params.side;
neig   = params.neig;
maxexp = params.maxexp;
minexp = params.minexp;

if(nargin<2)
  m = input('Please enter m=?');
end

if(nargin==0)
  u=input('Please enter u=?');
end

E = [];
m2     = m^2;
  
dims = 2.^(minexp:maxexp);
 
OPTS.disp = 0;
OPTS.maxit = 1000;
OPTS.issym = 1;
delta = 0;  
if(u>0)
  delta = log(exp(2*u)-1)/2;
end      

for n = dims,

  OPTS.p = min(lancz,n);	   % lancz < dim

  x = linspace(-side/2+delta,side/2+delta,n)';
  F=1+exp(4*u)+exp(2*x);
  G=1+exp(2*x);
  W = 1/3*exp(2*(u+x))./F./G*cosh(2*u);
  V = -1/3*exp(2*(u+x))./F./G.*(-3*exp(2*u)./G-6*cosh(2*u)*G./F)+m2*((1+tanh(x))/2).^2;
  
  [A,B] = AB_Matrices(x,V,W);  % build the sparse matrix 
  [PSI,Ecurr] =  eigs(A,B,neig,'SM',OPTS); % and diagonalize it with 'eigs'
  
  [Ecurr,ind]=sort(diag(Ecurr));
  PSI=PSI(:,ind);
  
  E = [E, Ecurr];           % save for Richardson extrapolation
end

Er = rich(E,order);
if params.plot=='y'
  subplot(211)
  fill(x-delta,V/6,'y','EraseMode','none');
  hold on
  plot(x-delta,W,'b','LineWidth',2);
  title('V/6, W')
  xlabel('\xi')
  legend('W')
  subplot(212)
  plot(x-delta,PSI.^2)
  title('\psi(\xi)')
  xlabel('\xi')
  
  for j=1:neig
    s{j}=['\psi_',num2str(j-1)];
  end
  
  legend(s)
    
  drawnow
  hold on
end

%-----------
function [A,B] = AB_Matrices(x,V,W)

side = max(x)-min(x);
n = length(x);
dx = side/(n-1);
e = ones(1,n)';
D2 = spdiags([-e e],0:1,n,n);
D2 = D2+D2';
D2(1,1) = -1;     % uncomment for Neumann b.b
D2(n,n) = -1;     % uncomment for Neumann b.b

D2 = D2/dx^2;     % the simplest (d/dx)^2 approx
                  % Richardson will  take care of
                  % finite dx errors
					
A = -D2+spdiags(V,0,n,n);		% factor 1/4 to be justified!
B =     spdiags(W,0,n,n);		% but it works!


%-------------
function y = rich(x,k)
% Ritorna il valore finale dello schema di Richardson
% a partire da dati valutati su un array n=2^m...2^M
% Se i dati sono in forma di matrice x(i,j), l'estrapolazione
% e' fatta per ogni "i" separatamente.
% Usage:   y=rich(x,k)
%        k=1,2,... dipende dalla forma dello sviluppo
%        x = a + b/n^k + c/n^2k + ...
% Ref. Davis  and Rabinowitz, "Numerical integration",
% Blaisdell Pub.Co., 1967

[n,m]=size(x);

if m==1,
  disp('Too few data, cannot extrapolate'); % need a table of values
  y = x;					% suitable for extrapolation
else
  fac = 2^k;
  for it = 1:m-1,
    for j = 1:m-it,	
      y(:,j) =  (fac * x(:,j+1) - x(:,j))/(fac-1);
    end
    x = y;
    fac = fac*2^k;
  end
end
y = y(:,1);
