function [la,resVec]=FEMEig(nodes,elem,edges,aFunc,bFunc,fFunc,eigVec,tol)
%[...] = FEMEig (...)
%  determine  eigenvalues la and eigenfunctions u for the BVP
%
%  div(a*grad u) - b*u = la*f*u   in domain
%                    u = 0        on Dirichlet section of the boundary
%              a*du/dn = 0        on Neumann section of the boundary 
%
%  la      = FEMEig(nodes,elem,edges,aFunc,bFunc,fFunc,eigVec,tol)
%  [la,ev] = FEMEig(nodes,elem,edges,aFunc,bFunc,fFunc,eigVec,tol)
%
% nodes elem, edges contains information about the mesh
%         see ReadMesh() for the description of the format
% aFunc bFunc fFunc  function files for the coefficient functions
%         may also be given as vectors or scalar values
% eigVec  is the initial guess for the eigenvectors
%         the number of columns determines the number of eigenvalues
%         to be computed
%         if a number n is given, then n eigenvalues will be computed
% tol is the tolerance for the relative error of the eigenvalues
%         if not given tol = 1e-5 is used as default
%
% la  is the vector containing the eigenvalues
% ev  is the matrix with the eigenvectors as columns
%
%see also ReadMesh, ShowMesh, FEMEquation, ShowSolution, FEMValue ...

if ((nargin<7)|(nargin>8)) 
  help("FEMEig");
  usage("FEMEig(nodes,elem,edges,aFunc,bFunc,fFunc,eigVec,tol)");
endif

if (nargin==7) tol=1e-5;endif

[gMat,gVec,n2d]=FEMEquation(nodes,elem,edges,aFunc,bFunc,fFunc,0,0);

if(size(eigVec)==[1,1])
  nn=eigVec;
%  eigVec=zeros(size(gMat)(1),eigVec);
%  eigVec(1:nn,1:nn)=eye(nn);
  eigVec=rand(size(gMat)(1),eigVec);
endif

if (nargout==1)
  la=SBEig(gMat,gVec,eigVec,tol);
endif

if (nargout==2)
  [la,ug]=SBEig(gMat,gVec,eigVec,tol);
  n=length(n2d);
  m=size(eigVec)(2);

  resVec=zeros(n,m);
  for k=1:n
    if n2d(k)>0
      resVec(k,:) = ug(n2d(k),:);
    endif
  endfor
endif
endfunction
