differential equations - Getting strange "not a polynomial" error from NDSolve

After answering this question, I decided to change it a bit (had some free time), and turned it into a Gross-Pitaevskii equation. Suddenly, the code that used to work now gives the error message:

CoefficientArrays::poly: "1/2 (-ψ$6110-ψ$6111)+ψ (1-
1/Sqrt[1+x$6089^2+y$6090^2]+0.1 Abs[ψ]^2) is not a polynomial.

I tried changing the potential for the Harmonic Oscillator, but still no use.

My code:

ClearAll["Global`*"];
g = 0.1;
X = 5;
V[x_, y_] := 1/2 (x^2 + y^2);
ini[x_, y_] := Exp[-(x^2 + y^2)];

bc = ψ[x, y] == 
    ini[x, y] /. {{x -> X}, {x -> -X}, {y -> X}, {y -> -X}};
sys = {-1/
       2 (D[ψ[x, y], {x, 2}] + 
        D[ψ[x, y], {y, 2}]) + (V[x, y] + 
        g Abs[ψ[x, y]^2]) ψ[x, y] == Etr ψ[x, y], bc};
sol = ParametricNDSolveValue[
  sys, ψ, {x, -X, X}, {y, -X, X}, {Etr}];
Plot3D[Abs[sol[1][x, y]], {x, -X, X}, {y, -X, X}, PlotRange -> All]

If I set g=0, I get some good results. Any other value and I get the error message.

Can anyone explain what's going on? How can it be fixed?

Answer

As mentioned in the comments above, currently "FiniteElement" method can't handle nonlinear PDE. Fortunately you're solving the PDE on a regular domain, so let's make a step backward i.e. try solving the equation with finite difference method (FDM). I'll use pdetoae for the generation of difference equations:

g = 0.1;
X = 5;
V[x_, y_] := 1/2 (x^2 + y^2);
ini[x_, y_] := Exp[-(x^2 + y^2)];
bc = ψ[x, y] == ini[x, y] /. {{x -> X}, {x -> -X}, {y -> X}, {y -> -X}};
eq = -1/2 (D[ψ[x, y], {x, 2}] + D[ψ[x, y], {y, 2}]) + (V[x, y] + 

       g Abs[ψ[x, y]^2]) ψ[x, y] == Etr ψ[x, y];

domain = {-X, X};
points = 25;
grid = Array[# &, points, domain];
difforder = 4;

(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
ptoafunc = pdetoae[ψ[x, y], {grid, grid}, difforder];


del = #[[2 ;; -2]] &;

ae = del /@ del@ptoafunc[eq];
aebc = MapAt[del, ptoafunc[bc], {{1}, {2}}];

var = Outer[ψ, grid, grid];

sol = 
   FindRoot[{ae, aebc} /. Etr -> 1, {#, 1} & /@ Flatten@var]; // AbsoluteTiming

(* {0.365270, Null} *)

ListPlot3D[var /. sol, PlotRange -> All, DataRange -> X]

Notice setting proper initial guess in FindRoot can be troublesome, but once again, fortunately your equation seems not to be the case.