[N,T] = NewMesh(0);                      % Generate square domain
[N,T]=RefineMesh(N,T);                   % Refine mesh once
err = zeros(2,5);                        % Error in the computed eigs
lambdas = zeros(2,5);                    % Value of eigs
for i=1:5,
    [N,T]=RefineMesh(N,T);               % Refine by factor of two
    [K,M,ind]=FEeig(N,T);                % Generate FEM matrices
    [V,D] = eigs(K,M,2,'SM');            % Generalized eig problem
    lambdas(:,i) = sort(diag(D));        % Sort eigs
    err(:,i) = lambdas(:,i)-[2;5]*pi^2   % Subtract exact solutions
    ratio = err(:,2:end)./err(:,1:end-1) % Ratio of errors (for conv. rate)
end;