Bilinearization with Mathematica - where to start?

I assume that by "bilinearize" you mean linearization relative to u[t, x, y] and psi[t, x, y]. This can be achieved using Series. Note, Series relative to more than one variable does does first linearization relative to the first variable. Then, the result is linearized relative to the second variable, e.t.c. That means, we get cross terms, that need to be eliminated. E.g.

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]};

To linearize, we do:

ser = Series[eq, {u[t, x, y], 0, 1}, {psi[t, x, y], 0, 1}];

This is a SeriesData object. For further work we need to get rid of the O[..] terms:

ser= ser // Normal;

To eliminate the cross terms, we first need to expand the expression and then set the cross terms to zero:

Expand[ser] /. u[t, x, y] psi[t, x, y] -> 0

Hint.

Assuming you have an approximation $(u_k,psi_k)$ then you can proceed linearly as

$$ cases{ mathcal{D}_1[u_{k+1}]+6(u_k psi_k +u_k(psi_{k+1}-psi_k)+psi_k(u_{k+1}-u_k))+(x^2+y^2)u_{k+1}=i mathcal{D}_2[u_{k+1}] -frac{i}{sqrt{2}}mathcal{D}_3[u_{k+1}]+psi_{k+1}=10^{-4}mathcal{D}_4[psi_{k+1}] } $$

solving for $(u_{k+1},psi_{k+1})$ to obtain the next approximation. Here $mathcal{D}_j[cdot]$ are differential operators.

Look what @Jens has done:

TaylorSeries[expr_, vars_?ListQ, start_?ListQ, order_?IntegerQ] := 
 Normal[Series[expr /. Thread[vars -> t*(vars - start) + start], {t, 0, order}]] /. t -> 1

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]};

TaylorSeries[eq, {u[t, x, y], psi[t, x, y]}, {0, 0}, 1]