function [K,M,ind]=FEeig(N,T);
% FEEIG generates the stiffness and mass matrices for the Laplace EVP
%   [K,M,IND]=FEeig(N,T); generates the stiffness matrix K and mass 
%   matrix M for the Laplace eigenvalue problem on a triangular mesh 
%   given by the list of triangles T and nodes N with Dirichlet boundary 
%   conditions for external edges. IND returns the node numbers of degrees
%   of freedom not on the Dirichlet boundary

n=size(T,1);
m=size(N,1);
bn=zeros(m,1);
ii=zeros(n,9);jj=zeros(n,9);                   % 9 entries for element matrice
Ka=zeros(n,9);Ma=zeros(n,9);  
[I,J]=ndgrid(1:3,1:3); I=I(:)'; J=J(:)';
for i=1:n,
  Ke=ComputeElementStiffnessMatrix([N(T(i,1),:); N(T(i,2),:); N(T(i,3),:)]);
  Me=ComputeElementMassMatrix([N(T(i,1),:); N(T(i,2),:); N(T(i,3),:)]);
  bn(T(i,1))=bn(T(i,1)) | T(i,4) | T(i,6);  % on the boundary
  bn(T(i,2))=bn(T(i,2)) | T(i,4) | T(i,5);
  bn(T(i,3))=bn(T(i,3)) | T(i,5) | T(i,6);
  ii(i,:)=T(i,I); jj(i,:)=T(i,J);             % assemble vectors
  Ka(i,:)=Ke(:)'; Ma(i,:)=Me(:)';             % for speed
end;
K=sparse(ii,jj,Ka);
M=sparse(ii,jj,Ma);

ind = find(bn==0);
K = K(ind,ind);
M = M(ind,ind);


% K=K';M=M';              % only for speed: deleting columns is faster!
% for i=1:m,
%   if bn(i)>0,
%     b(i)=feval(g,N(i,1),N(i,2));
%     K(:,i)=0;K(i,i)=1;
%     M(:,i)=K(:,i); M(i,:)=M(:,i)';
%   end;
% end;
% K=K';M=M'; 
% u=K\(M*b);