NDEigensystem to find eigenvalues and eigenfunctions of coupled differential equations:

Question

I would like to numerically solve the following system of coupled differential equations:

eq1=h/(2*ma) * PsiA''[x]+ VExt*PsiA[x]+ga*Abs[PsiA[x]]^2 + gab*Abs[PsiB[x]]^2*PsiA[x] == mua*PsiA[x]

eq2=h/(2*mb) * PsiB''[x]+ VExt*PsiB[x]+gb*Abs[PsiB[x]]^2 + gab*Abs[PsiA[x]]^2*PsiB[x] == mub*PsiB[x]

The following are known model parameters:

h ma mb ga gb gab

The term "VExt" is a known function of x, namely:

VExt=-P(Cos[3/2*x/L*2*Pi])^2

The boundary conditions are the periodic ones, i.e.

PsiA[x]==PsiA[x+L]

PsiB[x]==PsiB[x+L]

I would like to know how to find eigenvalues

mua mub

and eigenfunctions

PsiA[x]   PsiB[x]

Thanks in advance.

For example purposes, you can take:

h = 1; ma = 1; mb = 2; ga = 1; gab = 1; gb = 1; L = 4;

Answer 1

I will show the solution algorithm for some set of parameters.

h = 1; ma = 1; mb = 2; ga = 1; gab = 1; gb = 1; L = 4; VExt = \
1;
eq1 = h/(2*ma)*PsiA''[x] + VExt*PsiA[x] + ga*Abs[PsiA[x]]^2 + 
    gab*Abs[PsiB[x]]^2*PsiA[x] == mua*PsiA[x];

eq2 = h/(2*mb)*PsiB''[x] + VExt*PsiB[x] + gb*Abs[PsiB[x]]^2 + 
    gab*Abs[PsiA[x]]^2*PsiB[x] == mub*PsiB[x];

 bc1 = {PsiA[0] == 1, PsiA'[0] == 0};

bc2 = {PsiB[0] == 1, PsiB'[0] == 0};

 PSA = 
 ParametricNDSolveValue[{eq1, eq2, bc1, bc2}, 
  PsiA, {x, 0, L}, {mua, mub}];
PSB = ParametricNDSolveValue[{eq1, eq2, bc1, bc2}, 
  PsiB, {x, 0, L}, {mua, mub}];

FindRoot[{PSA[mua, mub][L] == 1, 
  PSB[mua, mub][L] == 1}, {{mua, -1}, {mub, -1}}]

Out[]= {mua -> -2.84339, mub -> -1.15353}

{Plot[PSA[-2.84339, -1.15353][x], {x, 0, L}, 
  AxesLabel -> {"x", "PsiA"}], 
 Plot[PSB[-2.84339, -1.15353][x], {x, 0, L}, 
  AxesLabel -> {"x", "PsiB"}]}

In the case of a periodic potential, it is possible to weaken the conditions at the boundary and use a parametric function with 4 parameters.

h = 1; ma = 1; mb = 2; ga = 1; gab = 1; gb = 1; L = 1; 
VExt[x_] := -(Cos[3*x/L*Pi])^2;
eq1 = h/(2*ma)*PsiA''[x] + VExt[x]*PsiA[x] + ga*Abs[PsiA[x]]^2 + 
    gab*Abs[PsiB[x]]^2*PsiA[x] == mua*PsiA[x];

eq2 = h/(2*mb)*PsiB''[x] + VExt[x]*PsiB[x] + gb*Abs[PsiB[x]]^2 + 
    gab*Abs[PsiA[x]]^2*PsiB[x] == mub*PsiB[x];
bc1 = {PsiA[0] == A0, PsiA'[0] == 0};

bc2 = {PsiB[0] == B0, PsiB'[0] == 0};
PSA = ParametricNDSolveValue[{eq1, eq2, bc1, bc2}, 
  PsiA, {x, 0, L}, {mua, mub, A0, B0}, 
  Method -> {"StiffnessSwitching", "NonstiffTest" -> False}, 
  MaxSteps -> Infinity]
PSB = ParametricNDSolveValue[{eq1, eq2, bc1, bc2}, 
  PsiB, {x, 0, L}, {mua, mub, A0, B0}, 
  Method -> {"StiffnessSwitching", "NonstiffTest" -> False}, 
  MaxSteps -> Infinity]

 FindRoot[{PSA[mua, mub, 1.9, 1.7][L] == 1.9, 
  PSB[mua, mub, 1.9, 1.7][L] == 1.7}, {{mua, -18}, {mub, -8}}, 
 Method -> "Secant"]

Out[]= {mua -> -17.8395, mub -> -7.76387}

 {Plot[PSA[-17.8395, -7.76387, 1.9, 1.7][x], {x, 0, L}, 
  AxesLabel -> {"x", "PsiA"}], 
 Plot[PSB[-17.8395, -7.76387, 1.9, 1.7][x], {x, 0, L}, 
  AxesLabel -> {"x", "PsiB"}]}