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

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_1ll 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_2gg 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.

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]

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

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