//  Eigen Values of Stelov Poincare Problem of Laplacian Operator
// in FreeFEM
int[int] l4=[1,1,1,1];
mesh Th=square(50,50,[x,y],label=l4);
fespace Vh(Th,P1);


macro Grad(u) [dx(u),dy(u)]//
real sigma = 1;//
varf vA(u,v)=int1d(Th)(-sigma*u*v)+int2d(Th)(Grad(v)'*Grad(u)) ; //'
varf vM(u,v)= int1d(Th)(u*v);
matrix A = vA(Vh,Vh,solver=sparsesolver);
matrix B = vM(Vh,Vh);
int nev = 20;
real[int] ev(nev);
Vh[int] eV(nev);
int k=EigenValue(A,B,sym=true,sigma=sigma,value=ev,vector=eV);
k=min(k,nev); //  some time the number of converged eigen value is  greater then nev..
cout << ev <<endl;
for(int i=0;i<k;++i)
 	plot(eV[i],cmm=" vp " + ev[i],wait=1,dim=3,fill=1);
