Matching the solutions of diff. equations from forward and backward in some point

Question

I am trying to solve two coupled non-linear differential equations for $F(r)$ and $h(r)$: $$ \begin{aligned} F''-F(F^2-1)/r^2- Fh^2&=0 \\ h''+2h'/r-2F^2h/r^2+\beta^2/2 h(1-h^2)&=0 \end{aligned} $$ I know the solutions' behaviour near the origin (left): $$ \begin{aligned} F&=1+a r^2+... \\ h&=b r+... \end{aligned} $$ and at the infinity (right): $$ \begin{aligned} F&=Ae^{-r}+... \\ h&=1-B e^{-\beta r}/r+... \end{aligned} $$ where $a$, $b$, $A$ and $B$ are free parameters, $\beta$ is fixed.

To solve the system of equations I do shooting forward from $r=r_1\ll 1$:

r1 = 0.01;
r2 = 10;
mid = 4;
beta=1;

profFun1 = 
 ParametricNDSolveValue[{F''[r] - F[r] (F[r]^2 - 1)/r^2 - 
     F[r] h[r]^2 == 0, 
   h''[r] + 2/r h'[r] - 2 F[r]^2 h[r]/r^2 + 
     beta^2/2 h[r] (1 - h[r]) (1 + h[r]) == 0, 
   F[r1] == 1 + a r1^2 + 1/10 (3 a^2 + b^2) r1^4, 
   F'[r1] == 2 a r1 + 4/10 (3 a^2 + b^2) r1^3, 
   h[r1] == b r1 + b/10 (4 a - beta) r1^3, 
   h'[r1] == b + 3 b/10 (4 a - beta) r1^2}, {F, h}, {r, r1, mid}, {a, 
   b}, Method -> "ExplicitRungeKutta"]

and backward from $r=r_2\gg 1$ using the substitution $u=1/r$:

profFun2 = 
 ParametricNDSolveValue[{ 
   F1''[u] + 2/u F1'[u] - F1[u] (F1[u]^2 - 1)/u^2 - 
     F1[u] h1[u]^2/u^4 == 0, 
   h1''[u] - 2 F1[u]^2 h1[u]/u^2 + 
     beta^2/2 /u^4 h1[u] (1 - h1[u]) (1 + h1[u]) == 0, 
   F1[r1] == A Exp[-1/r1], F1'[r1] == A Exp[-1/r1] 1/r1^2, 
   h1[r1] == 1 + B Exp[-beta /r1] r1, 
   h1'[r1] == B (Exp[-beta /r1] + beta Exp[-beta/r1]/r1)}, {F1, 
   h1}, {u, r1, 1/mid}, {A, B}]

So, I get the solutions (for the left and the right):

FFl[a_, b_, r_] := profFun1[a, b][[1]][r]
hhl[a_, b_, r_] := profFun1[a, b][[2]][r]
FFr[A_, B_, r_] := profFun2[A, B][[1]][1/r]
hhr[A_, B_, r_] := profFun2[A, B][[2]][1/r]

Now I need to find such a, b, A, B that in a point r=mid the functions from the left and the right are matched smoothly. For this I also create derivatives of the functions in the point r=mid

FFld[a_, b_] := D[FFl[a, b, r], r] /. r -> mid
hhld[a_, b_] := D[hhl[a, b, r], r] /. r -> mid
FFrd[A_, B_] := D[FFr[A, B, r], r] /. r -> mid
hhrd[A_, B_] := D[hhr[A, B, r], r] /. r -> mid

For matching I have these four equations (for functions and their derivatives)

{FFl[a, b, mid] == FFr[A, B, mid], hhl[a, b, mid] == hhr[A, B, mid], 
FFld[a, b] == FFrd[A, B], hhld[a, b] == hhrd[A, B]}

So, my target is to find these four parameters to match the function from the left and the right. For my further calculations, I need just functions $F(r)$ and $h(r)$ independent on the free parameters.

I have tried something like this

FindRoot[{FFl[a, b, mid] == FFr[A, B, mid], hhl[a, b, mid] == hhr[A, B, mid], 
FFld[a, b] == FFrd[A, B], hhld[a, b] == hhrd[A, B]}, {{a,-0.3}, {b,0.6},{A,1.},{B,0.05}}]

but it does not work. However, I had already done this algorithm with FindRoot for one non-linear differential equation (with two free parameters) and it worked.

I would really appreciate suggestions, how I can get these parameters. I also used Solve instead of FindRoot, it does not work too. Maybe I should "say" something additional to FindRoot for helping it. Maybe there are some better alternatives to FindRoot for such task?

I would be also very thankful for any other ideas for solving such a system of differential equations. I recently got another system with four equations and accordingly with eight free parameters. I found this idea to solve such differential equations by matching the solutions from the left and the right side in some article, but do not know how to realise it in Mathematica. It would be great for me to learn it.

Thank you very much for reading.

Answer 1

I would be also very thankful for any other ideas for solving such a system of differential equations.

Then why not new-in-v12 nonlinear FEM of NDSolve?:

r1 = 0;
r2 = 9;
beta = 1;

sol=NDSolveValue[{F''[r] - F[r] (F[r]^2 - 1)/r^2 - F[r] h[r]^2 == 0, 
     h''[r] + 2/r h'[r] - 2 F[r]^2 h[r]/r^2 + 
         beta^2/2 h[r] (1 - h[r]) (1 + h[r]) == 0, 
     F[r1] == 1, h[r1] == 0, F[r2] == 0, h[r2] == 1}, {F, h}, {r, r1, r2}, 
 Method -> FiniteElement, InitialSeeding -> {F[r]==1-r/r2,h[r]==r/r2}]

Plot[sol[t] // Through // Evaluate, {t, r1, r2}, PlotRange -> All]

Notice I've made r2 a bit smaller. With better initial guess (it's set by the InitialSeeding option) one should be able to set larger r2.

"But what's wrong with my original attempt?" There're several issues:

1. Evaluation order isn't controlled properly. Just execute FFl[a, b, mid] == FFr[A, B, mid] and observe the output, you'll see what's wrong. Add ?NumericQ to proper positions or add Evaluated -> False to FindRoot will resolve the problem.

2. Value of mid isn't proper, it should be set to 1, as suggested by the asymptotic solutions.

3. You've used r1 rather than r2 in profFun2. This is equivalent to setting r2=100, which is way too large.

After correcting all these, one can obtain:

{a -> -0.340902, b -> 0.73182, A -> 3.17522, B -> -1.90734}

Aside from these critical problems, your implementation is unnecessarily verbose. The following is mine, compare it with yours carefully:

r1 = 1/100;
r2 = 10;
mid = 1;
β = 1;

eq = {F''[r] - F[r] (F[r]^2 - 1)/r^2 - F[r] h[r]^2 == 0, 
      h''[r] + 2/r h'[r] - 2 F[r]^2 h[r]/r^2 + β^2/2 h[r] (1 - h[r]) (1 + h[r]) == 0};

FL[r_] = 1 + a r^2 + 1/10 (3 a^2 + b^2) r1^4;
hL[r_] = b r + b/10 (4 a - β) r^3;
FR[r_] = A Exp[-r];
hR[r_] = 1 + B Exp[-β r] /r;

sys1 = {eq, F[r1] == FL[r1], F'[r1] == FL'[r1], h[r1] == hL[r1], h'[r1] == hL'[r1]};
sys2 = {eq, F[r2] == FR[r2], F'[r2] == FR'[r2], h[r2] == hR[r2], h'[r2] == hR'[r2]};

varmid = {F, h, F', h'}[mid] // Through

profFun1 = ParametricNDSolveValue[sys1, varmid, {r, r1, mid}, {a, b}];
profFun2 = ParametricNDSolveValue[sys2, varmid, {r, mid, r2}, {A, B}];

solrule = FindRoot[
  profFun1[a, b] - profFun2[A, B], {{a, -0.003}, {b, 0.738}, {A, 1}, {B, 0.05}}]

ListLinePlot@
 Flatten@{NDSolveValue[sys1 /. solrule, {F, h}, {r, r1, mid}], 
          NDSolveValue[sys2 /. solrule, {F, h}, {r, mid, r2}]}