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