Ok I'm sorry I was wrong... With Joachims comment I think it's actually quite easy...
Code:
from netgen.geom2d import unit_square
from ngsolve import *
mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))
h1=H1(mesh,order=5,dirichlet="bottom|left|right|top",complex=True)
u,v = h1.TnT()
a = BilinearForm(h1)
a += SymbolicBFI(grad(u)*grad(v))
b = BilinearForm(h1)
b += SymbolicBFI(u*v)
a.Assemble()
b.Assemble()
evecs = GridFunction(h1,multidim=50)
vlam = ArnoldiSolver(a.mat,b.mat,h1.FreeDofs(),evecs.vecs,1)
Draw(evecs)
# get a non multidim GF out of the multidim one
evec0 = GridFunction(h1)
evec0.vec.data = evecs.vecs[0]
number = NumberSpace(mesh,complex=True)
fes = FESpace([h1,number])
(u,unum), (v,vnum) = fes.TnT()
# The orthogonality condition is satisfied because vnum = 1 means that evec0 * u = 0
# unum incorporates the evec0 part of the right hand side, if the argument of f is orthogonal to
# evec0, then unum should be 0
alam = BilinearForm(fes)
alam += SymbolicBFI(grad(u) * grad(v) - vlam[0] * u * v + unum * v * evec0)
alam += SymbolicBFI(evec0 * u * vnum)
f = LinearForm(fes)
f += SymbolicLFI(1*v)
alam.Assemble()
f.Assemble()
sol = GridFunction(fes,"solution")
sol.vec.data = alam.mat.Inverse(fes.FreeDofs()) * f.vec
Draw(sol.components[0])
