tol = 1e-8;
ngrid = 20;


G = numgrid('L',ngrid+1);
h = 2/ngrid;
A = delsq(G)/h^2;
n = size(A,1);

lower = sort(eigs(A,2,'SM'));
upper = sort(eigs(A,2,'LM'));

lambdamin = lower(1);
lambda2 = lower(2);
lambdamax = upper(2);
lambdapen = upper(1);

maxiter = 40;
q0 = rand(n,1);
alpha = zeros(maxiter,1);
beta = zeros(maxiter-1,1);
q0 = q0/norm(q0);
qlast = zeros(n,1);

ritz = zeros(maxiter);

for i=1:maxiter,
  r = A*q0;
  alpha(i) = q0'*r;
  r = r - alpha(i)*q0;
  if (i>1),
     r = r - beta(i-1)*qlast;
  end;
  beta(i) = norm(r);
  if (beta(i) < tol),
    break;
  end;
  qlast = q0;
  q0 = r/beta(i);

  if (i==1),
    ritz(1,1) = alpha(1);
  else
    T = diag(alpha(1:i),0) + diag(beta(1:i-1),1) + diag(beta(1:i-1),-1);
    ritz(1:i,i) = sort(eig(T));
  end;
end;

x1 = (lambdamax - lambda2)/(lambda2 + lambdamax - 2*lambdamin);
x2 = (lambdapen - lambdamin)/(2*lambdamax - lambdapen - lambdamin);
bound1 = acosh(1/x1);
bound2 = acosh(1/x2);

%semilogy(1:i,ritz(1,1:i)-lambdamin,'kv-',1:i,1./(cosh((1:i)*bound1)).^2, 'k--',...
%         1:i,lambdamax-diag(ritz(1:i,1:i)),'b^-',1:i,1./(cosh((1:i)*bound2)).^2,'b--');
semilogy(1:i,ritz(1,1:i)-lambdamin,'kv-',...
         1:i,lambdamax-diag(ritz(1:i,1:i)),'b^-',...
         1:i,(lambdamax-ritz(1,1))./(cosh((1:i)*bound1)).^2, 'k--');

legend('Smallest','Largest','Bound');
