Bilinearization with Mathematica – where to start?

I tried to bilinearize the two equations:

eq = {Laplacian(psi(t, x, y), {x, y}) + 
     6 u(t, x, y) psi(t, x, y) + (x^2 + y^2) u(t, x, y) == 
    I D(u(t, x, y), 
      t), -(I /Sqrt(2) (D(u(t, x, y), x) + D(u(t, x, y), y)) + 
       psi(t, x, y)) == 1/10000 D(psi(t, x, y), t)};

and I have looked at some posts here on stackexchange, but found no posts that considered the matter of bilinearization. Does anyone have a tip of where to start?

nssm2 = NonlinearStateSpaceModel(eq, psi(t,x,y), u(t,x,y), t)

but to no avail.
I read also this:

https://library.wolfram.com/infocenter/Articles/2831/

But I couldn’t find a hint there either.

Thanks