# NDEigensystem to find eigenvalues and eigenfunctions of coupled differential equations:

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]


For example purposes, you can take:

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

• Indicate the particular case of functions and parameters. I will show the solution algorithm. Oct 23, 2018 at 15:13
• I've changed my post specifying the form of VExt. Concerning model parameters, what you've chosen is very good. Can you please modify your answer accordingly? Oct 23, 2018 at 15:39
• Can you also briefly explain and comment your code. In particular, I don't understand the role of bc1 and bc2. They seem to be further constraints which I have no reasons to introduce... Oct 23, 2018 at 15:41
• These are the usual conditions for such tasks. To find the eigenvalues, three boundary conditions must be specified, not two. Oct 23, 2018 at 16:11
• Periodic solutions may not exist for every potential. Oct 23, 2018 at 16:45

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