As initial step to compute ground state we can use LDG method as follows
Clear["Global`*"]
Get["NumericalDifferentialEquationAnalysis`"];
L = 10; \[Beta] = 100; \[Omega] = 1; muTF =
1/2 (3 \[Omega] \[Beta]/2)^(2/3);
f = {2 (- mu p[t] + 1/2*\[Omega]^2 t^2 p[t] + \[Beta] p[t]^2 p[t]),
v[t]}; ini = {0, b1}; bc1 = 0;
LDGODEs[M0_, nn_, ns_, np_, f_, ini_, tmax_, tmin_] :=
Module[{dx = (tmax - tmin)/nn, A = Array[a, {M0 + 1, nn, 2}]},
xl = Table[tmin + l*dx, {l, 0, nn}]; UT[m_, t_] := EulerE[m, t];
psi[m_, k_, t_] :=
Piecewise[{{UT[
m, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])],
xl[[k]] <= t <= xl[[k + 1]]}, {0, True}}];
g = Table[
GaussianQuadratureWeights[np, xl[[i]], xl[[i + 1]]], {i, nn}];
dp = Table[
D[UT[m, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])],
t], {k, nn}];
rul = {p[t] ->
Sum[a[m, k, 2] psi[m - 1, k, t], {m, 1, M0 + 1}, {k, 1, nn}],
v[t] -> Sum[
a[m, k, 1] psi[m - 1, k, t], {m, 1, M0 + 1}, {k, 1, nn}]};
eq = Flatten[
Table[-Sum[
a[i + 1, n,
ks] (Table[(psi[i, n, t] If[j == 0, 0,
dp[[n]] /. m -> j]), {t, g[[n]][[All, 1]]}] .
g[[n]][[All, 2]]), {i, 0,
M0}] - (Table[(f[[ks]] /. rul) psi[j, n, t], {t,
g[[n]][[All, 1]]}] . g[[n]][[All, 2]]), {j, 0, M0}, {n, 1,
nn}, {ks, 1, 2}] +
Table[Sum[
a[i + 1, n,
ks] (psi[i, n, xl[[n + 1]]] psi[j, n, xl[[n + 1]]]), {i, 0,
M0}] - Sum[
If[n < 2, ini[[ks]]/(M0 + 1),
a[i + 1, n - 1, ks] psi[i, n - 1, xl[[n]]]] psi[j, n,
xl[[n]]], {i, 0, M0}], {j, 0, M0}, {n, 1, nn}, {ks,
2 ns}]];
eq1 = (p[t] - bc1) /. rul /. t -> tmax;
eqn = Join[eq, {eq1}];
var = Join[Flatten[A], {b1, mu}];
energy =
Sum[Table[((1/2 v[t]^2 +
1/2 \[Omega]^2 t^2 p[t]^2 + \[Beta]/2 p[t]^4) /. rul), {t,
g[[n]][[All, 1]]}] . g[[n]][[All, 2]], {n, 1, nn}];
norm = Sum[
Table[(( p[t]^2) /. rul), {t, g[[n]][[All, 1]]}] .
g[[n]][[All, 2]], {n, 1, nn}] == 1/2;
con = Join[{norm}, Table[eqn[[i]] == 0, {i, Length[eqn]}]];
sol1 = NMinimize[{energy, con}, var]; sol1]
Solution
M0 = 2; nn = 10; ns = 1; np = 5; tmax = L; tmin = 0;
ldgsol =
LDGODEs[M0, nn, ns, np, f, ini, tmax, tmin]; // AbsoluteTiming
Please note, that we minimize energy to compute $\mu _g$, Finally we have
mu /. sol1[[2]]
Out[]= 14.1341
This value we can compare to theoretical value muTF = 14.1155
. Wave function we should extend on $-L\le x\le L$ as
lst = Table[{t,
Evaluate[
Sum[a[i + 1, k, 2] psi[i, k, Abs[t]], {i, 0, M0}, {k, nn}] /.
sol1[[2]]]}, {t, -L, L, .0101}];
\[Phi]g = Interpolation[Join[lst, {{10, 0}}], InterpolationOrder -> 4];
Compare this function with theoretical one (red line)
Show[Plot[\[Phi]g[x], {x, -L, L}],
Plot[If[muTF > x^2/2, Sqrt[(muTF - x^2/2)/100], 0], {x, -10, 10},
PlotStyle -> {Red, Dashed}]]

Using LDG method we can compute several excitations, but problem is correct on $-\infty <x<\infty$ only, while on limited interval we have some mixture of states. For example
ns = 2;μg=mu/.sol1[[2]]; x1[t_] = Table[Symbol["x1" <> ToString[i]][t], {i, 1, ns}];
v1[t_] =
Table[Symbol["v1" <> ToString[i]][t], {i, 1, ns}]; bc = {0, 0}; f1 =
Join[{2 (0.5 \[Omega]^2 x^2 vl[x] +
2 \[Beta] Abs[\[Phi]g[x]]^2 vl[x] - \[Mu]g vl[
x] - \[Beta] \[Phi]g[x]^2 ul[x] + mu1 vl[x]),
2 (0.5 \[Omega]^2 x^2 ul[x] +
2 \[Beta] Abs[\[Phi]g[x]]^2 ul[x] - \[Mu]g ul[
x] - \[Beta] \[Phi]g[x]^2 vl[x] - mu1 ul[x])} /. {ul[x] ->
x11[t], vl[x] -> x12[t], \[Phi]g[x] ->
If[Abs[t] <= L, \[Phi]g[t], 0]} /. x -> t, v1[t]];
LDGODE[M0_, nn_, ns_, np_, f_, ini_, tmax_, tmin_] :=
Module[{dx = (tmax - tmin)/nn, A = Array[a, {M0 + 1, nn, 2 ns}]},
xl = Table[tmin + l*dx, {l, 0, nn}];
UT[m_, t_] := BernoulliB[m, t];
psi[m_, k_, t_] :=
Piecewise[{{UT[
m, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])],
xl[[k]] <= t <= xl[[k + 1]]}, {0, True}}];
g = Table[
GaussianQuadratureWeights[np, xl[[i]], xl[[i + 1]]], {i, nn}];
dp = Table[
D[UT[m, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])],
t], {k, nn}];
rul1 =
Join[Table[
x1[t][[i]] ->
Sum[a[m, k, i + ns] psi[m - 1, k, t], {m, 1, M0 + 1}, {k, 1,
nn}], {i, ns}],
Table[v1[t][[i]] ->
Sum[a[m, k, i] psi[m - 1, k, t], {m, 1, M0 + 1}, {k, 1,
nn}], {i, 1, ns}]];
eq = Flatten[
Table[-Sum[
a[i + 1, n,
ks] (Table[(psi[i, n, t] If[j == 0, 0,
dp[[n]] /. m -> j]), {t, g[[n]][[All, 1]]}] .
g[[n]][[All, 2]]), {i, 0,
M0}] - (Table[(f[[ks]] /. rul1) psi[j, n, t], {t,
g[[n]][[All, 1]]}] . g[[n]][[All, 2]]), {j, 0, M0}, {n, 1,
nn}, {ks, 2 ns}] +
Table[Sum[
a[i + 1, n,
ks] (psi[i, n, xl[[n + 1]]] psi[j, n, xl[[n + 1]]]), {i, 0,
M0}] - Sum[
If[n < 2, ini[[ks]]/(M0 + 1),
a[i + 1, n - 1, ks] psi[i, n - 1, xl[[n]]]] psi[j, n,
xl[[n]]], {i, 0, M0}], {j, 0, M0}, {n, 1, nn}, {ks,
2 ns}]];
eqb = {(x1[t][[1]] - bc[[1]]), x1[t][[2]] - bc[[2]]} /. rul1 /.
t -> tmax;
var = Join[Flatten[A], {b1, b2, mu1}];
norm1 = {Sum[
Table[(( x11[t]^2 - x12[t]^2) /. rul1), {t,
g[[n]][[All, 1]]}] . g[[n]][[All, 2]], {n, 1, nn}] == 1/2};
eqn = Join[eq, eqb];
sol2 =
FindRoot[Join[Table[eqn[[i]] == 0, {i, Length[eqn]}], norm1],
Table[{var[[i]], 1/10}, {i, Length[var]}]]; sol2];
Solution
ini = Join[{b1, b2}, {0, 0}]; M0 = 4; nn = 40; ns = 2; np = 5;
tmax = 9.9863; tmin = 0;
ldgsol2 =
LDGODE[M0, nn, ns, np, f1, ini, tmax, tmin]; // AbsoluteTiming
It should be first eigenvalue $\omega_1=1$, but we have instead
mu1 /. sol2
Out[]= 1.00091
Nevertheless, we have norm 2 norm1[[1, 1]] /. sol2
=1.
Theoretical solution for eigenfunction is given by
u0[x_] := 1/Sqrt[2] (x \[Phi]g[x] - \[Phi]g'[x]);
v0[x_] := 1/Sqrt[2] (x \[Phi]g[x] + \[Phi]g'[x])
Numerical solution is
pl1 = Table[
Plot[Evaluate[
Sum[a[i + 1, k, s] psi[i, k, t], {i, 0, M0}, {k, nn}] /.
sol2], {t, tmin, tmax}, PlotRange -> All, PlotStyle -> {Red},
Exclusions -> None], {s, ns + 1, 2 ns}]
Compare to theory they look like
pl = {Plot[u0[x], {x, 0, Sqrt[2 muTF] - .07}],
Plot[v0[x], {x, 0, Sqrt[2 muTF] - .07}]}
Show[pl1, pl]

{gvals, gfuns} = NDEigensystem[...
is that the expression is not linear in\[Phi]
. Also, stick to default methods unless you know what you are doing. $\endgroup$f[_,_]:=Abort[]; NDEigensystem[f''[x],f[x],{x,-1,1},1]
which leads to an abort in V12.3, which I find surprising (anyone?). For this reason, I suggest you drop the\[Phi][n_,x_]:=...
definition in your code, do not forget to clear old definitions or restart your kernel, and perhaps use another name for it that does not interfere with the function\[Phi]
in yourNDEigensystem
. $\endgroup$NDEigensystem
since problem is nonlinear any way due to constraint $\int {(u^2-v^2)dx=1}$. I think that we can try to solve nonlinear problem with usingNMinimize
to optimize energy. $\endgroup$