% Chebfun code to compute the eigenvalues of the Laplacian on L^2[0,pi]
% with homogeneous Dirichlet boundary conditions by applying
% the Arnoldi (Lanczos) algorithm directly to the inverse operator.

% Rice Eigenvalue Clinic, August 2009

 m = 10;                         % maximum number of iterations
 v0 = chebfun('x.^0',[0 pi]);    % starting vector - bad: orthogonal to even eigenfuns
 v0 = chebfun('exp(x)',[0 pi]);  % starting vector
 xfun = chebfun('x',[0 pi]);     % the function "x"
 V = v0/norm(v0);                % normalize the starting vector
 H = zeros(m+1,m);               % H will contain the upper Hessenberg matrix

 for k=1:m                       % Arnoldi process
    w = -cumsum(cumsum(V(:,k))); %            apply the inverse operator
    w = w - (w(0) + (w(pi)-w(0))*xfun)/pi;  % and patch up boundary conditions
    for j=1:k                    %    Gram-Schmidt
        H(j,k) = V(:,j)'*w;
        w = w-H(j,k)*V(:,j);
    end
    H(k+1,k) = norm(w);
    V(:,k+1) = w/H(k+1,k);       %    append new basis vector to V

    fprintf('\n-----------------------------------------\n')
    fprintf('step %d\n', k);
    fprintf('-----------------------------------------\n')
    fprintf('      Ritz value       relative error \n')
    fprintf('-----------------------------------------\n')
    ew = sort(1./eig(H(1:k,1:k)));
    for j=1:k
       fprintf('%20.15f    %12.5e\n', ew(j), (j^2-ew(j))/(j^2))
    end
 end

 fprintf('\nDeparture from orthogonality of Arnoldi basis functions: %12.5e\n\n', ...
     norm(V'*V-eye(m+1)))