if you use the Arnoldi solver of NGSolve then you store the modes in the multidim gridfunction:
Code:
from netgen.geom2d import unit_square
from ngsolve import *
mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))
fes = H1(mesh,order=5,dirichlet="bottom|left|right|top",complex=True)
u,v = fes.TrialFunction(), fes.TestFunction()
a = BilinearForm(fes)
a += SymbolicBFI(grad(u)*grad(v))
b = BilinearForm(fes)
b += SymbolicBFI(u*v)
a.Assemble()
b.Assemble()
u = GridFunction(fes,multidim=50)
vlam = ArnoldiSolver(a.mat,b.mat,fes.FreeDofs(),u.vecs,1)
Draw(u)
Then you can visualize the different eigenfunctions with Visual -> and change multidim component
