Boundary value problem, multiple dimensional shooting, coupled eigenvalue problem

Question

Following the one dimensional boundary value problem here, I would like to understand the easiest way to solve a BVP for a coupled system. In the 1D case, BVP can be converted to an initial value problem, and solved using shooting method. It was easily automated using bisection method.

On the other hand, when there are more than one boundary value to guess, bisection is simply not possible. Therefore, the question is, what would be the most efficient method for such problems.

To show the problem concretely, I have the following ODE system.

\begin{align} ( \frac{\tilde \mu_1^2}{B} -1 ) \tilde \Phi_1^2 + \frac{1}{A} {\tilde \Phi_1^{\prime 2}} + \frac{1}{2}\tilde \lambda_1 \tilde \Phi_1^4 + \frac{1}{2}\tilde \lambda_{12} \tilde \Phi_1^2 \tilde \Phi_2^2 +\frac{\tilde \mu_2^2}{B} \tilde \Phi_2^2 + \frac{1}{A} {\tilde \Phi_2^{\prime 2}} - \frac{A'}{\tilde r A^2} + \frac{1}{\tilde r^2 A} - \frac{1}{\tilde r^2} &= 0, \cr ( \frac{\tilde \mu_1 ^2 }{B} + 1 )\tilde \Phi_1^2 + \frac{1}{A} \tilde \Phi_1^{\prime 2} - \frac{1}{2}\tilde \lambda_1 \tilde \Phi_1^4 - \frac{1}{2}\tilde \lambda_{12} \tilde \Phi_1^2 \tilde \Phi_2^2 + \frac{\tilde \mu_2 ^2 }{B} \tilde \Phi_2^2 + \frac{1}{A} \tilde \Phi_2^{\prime 2} -\frac{B'}{\tilde r A B} - \frac{1}{\tilde r^2 A} + \frac{1}{\tilde r^2} &= 0, \cr \frac{1}{A}{\tilde \Phi_1^{\prime\prime }} +\left (\frac{\tilde \mu_1^2}{B} + 1 \right )\tilde \Phi_1 + \tilde \Phi'_1 \left (\frac{B'}{2 AB} - \frac{A'}{2A^2} + \frac{2}{A \tilde r}\right ) - \tilde \lambda_1 \tilde \Phi_1^3 - \frac{\tilde \lambda_{12}}{2} \tilde \Phi_2^2 \tilde \Phi_1 &= 0, \cr \frac{1}{A}{\tilde \Phi_2^{\prime\prime }} +\left (\frac{\tilde \mu_2^2}{B} \right )\tilde \Phi_2 + \tilde \Phi'_2 \left (\frac{B'}{2 AB} - \frac{A'}{2A^2} + \frac{2}{A \tilde r}\right ) - \frac{\tilde \lambda_{12}}{2} \tilde \Phi_1^2 \tilde \Phi_2 &= 0, \end{align} where $\tilde \Phi_1(\tilde r), \tilde \Phi_2(\tilde r), A(\tilde r), B(\tilde r)$ are functions to be solved, with boundary condition \begin{align} A(\infty) & = 1, \cr B(\infty) & = 1, \cr \tilde \Phi_1'(0) & = 0, \cr \tilde \Phi_1(\infty) & = 0, \cr \tilde \Phi_2'(0) & = 0, \cr \tilde \Phi_2(\infty) & = 0, \end{align} So there are $four$ functions, $six$ undetermined coefficients related to all the derivatives. On the other hand, there are $four$ equations constraints, and $six$ boundary values, so it should be solvable. There should be a class of $(\tilde \mu, \tilde \Phi)$, among which I want to select the 'ground state', $i.e.$ those $\tilde \Phi$'s that do not have zeros at finite $\tilde r$.

However, FEM does not handle it as it is not a linear coupled ODE system. Using NDSolve directly takes too long and never spits out any results, although there is no error thrown out either. Any idea?

Code:

Clear[GRHunterPhi4twoS2]
GRHunterPhi4twoS2[Omega_, Omega2_, \[CapitalLambda]1_, \[CapitalLambda]12_, xEnd1_] := 
Module[{}, Arule = {A[x] -> (1 - 2*(M[x]/x))^(-1), Derivative[1][A][x] -> 
   D[(1 - 2*(M[x]/x))^(-1), x]}; 
eq1 = -Derivative[1][M][x] + x^2*((1/2)*(Omega^2/B[x] - 1)*\[Sigma][x]^2 + (1/2)*(Omega2^2/B[x])*
       \[Sigma]2[x]^2 + (\[CapitalLambda]1/4)*\[Sigma][x]^4 + (\[CapitalLambda]12/4)*\[Sigma][x]^2*\[Sigma]2[x]^2 + Derivative[1][\[Sigma]][x]^2/2/A[x] + 
      Derivative[1][\[Sigma]2][x]^2/2/A[x]) == 0 /. Arule; 
eq2 = Derivative[1][B][x]/A[x]/B[x]/x - (1/x^2)*(1 - 1/A[x]) - (Omega^2/B[x] + 1)*\[Sigma][x]^2 - 
    (Omega2^2/B[x])*\[Sigma]2[x]^2 + (\[CapitalLambda]1/2)*\[Sigma][x]^4 + (\[CapitalLambda]12/2)*\[Sigma][x]^2*\[Sigma]2[x]^2 - 
    Derivative[1][\[Sigma]][x]^2/A[x] - Derivative[1][\[Sigma]2][x]^2/A[x] == 0 /. Arule; 
eq3 = Derivative[2][\[Sigma]][x] + (2/x + Derivative[1][B][x]/2/B[x] - Derivative[1][A][x]/2/A[x])*
     Derivative[1][\[Sigma]][x] + A[x]*((Omega^2/B[x] + 1)*\[Sigma][x] - \[CapitalLambda]1*\[Sigma][x]^3 - 
      (\[CapitalLambda]12/2)*\[Sigma]2[x]^2*\[Sigma][x]) == 0 /. Arule; 
eq4 = Derivative[2][\[Sigma]2][x] + (2/x + Derivative[1][B][x]/2/B[x] - Derivative[1][A][x]/2/A[x])*
     Derivative[1][\[Sigma]2][x] + A[x]*((Omega2^2/B[x])*\[Sigma]2[x]) - (\[CapitalLambda]12/2)*\[Sigma]2[x]*\[Sigma][x]^2 == 0 /. 
  Arule; bc = {M[xStart1] == epsilon, Derivative[1][\[Sigma]][xStart1] == epsilon, 
  Derivative[1][\[Sigma]2][xStart1] == epsilon, B[xEnd1] == 1, \[Sigma][xEnd1] == epsilon, 
  \[Sigma]2[xEnd1] == epsilon}; sollst1 = NDSolveValue[Flatten[{eq1, eq2, eq3, eq4, bc}], 
  {M, B, \[Sigma], \[Sigma]2}, {x, xStart1, xEnd1}, Method -> "StiffnessSwitching", 
  WorkingPrecision -> WorkingprecisionVar, AccuracyGoal -> accuracyVar, 
  PrecisionGoal -> precisionGoalVar, MaxSteps -> Infinity]]
precisionVar = 40; 
psiEndTolerance = 10^(-5); 
epsilon = 10^(-6); 
xStart1 = 10^(-5); 
xEnd1 = 30; 
WorkingprecisionVar = 23; 
accuracyVar = 9; 
precisionGoalVar = 9; 
Off[NDSolveValue::precw]
tmpSol = GRHunterPhi4twoS2[2, 2, 100, 50, xEnd1]
Plot[tmpSol[[4]][tmpx], {tmpx, xStart1, xEnd1}, PlotRange -> All]

P.S: my thought is somehow to convert the above equation set into a matrix representation. Then, maybe one can make use of the built-in Eigensystem[]. In this sense, with the boundary condition fulfilled, it can also be treated as an eigenvalue problem for $(\tilde \mu_1, \tilde \Phi_1)$ and $(\tilde \mu_2, \tilde \Phi_2)$. However, the non-linear terms $\tilde \Phi^n$ make this hard as well.