How to improve the accuracy of FindRoot or NDSolveValue at extreme conditions in my equation system?

Question

This is a follow up question to my earlier one. I have posted a question related to a DAE system. Here is my earlier question and the equation system now I am asking for an answer is

 L = 10; l = L*(3/1000); DD = 0.01; mu = 1;
     sol = 
         ParametricNDSolve[{ϕ1''[s] - (T - F)*Sin[ϕ1[s]] == 
            0, ϕ2''[s] - T*Sin[ϕ2[s]] == 0, ϕ1[0] == 
            0, ϕ2'[L] == 0, ϕ1[l] == ϕ2[l] == 
            ArcTan[1/mu]}, {ϕ1, ϕ2}, {s, 0, L}, {T, F}]; 
        m1[s_, T_, F_] := 
         Evaluate[Evaluate[ϕ1[T, F] /. sol][s]]; 
        m2[s_, T_, F_] := 
         Evaluate[Evaluate[ϕ2[T, F] /. sol][s]]; 
        dm1[s_, T_, F_] := Evaluate[D[m1[s, T, F], s]]; 
        dm2[s_, T_, F_] := Evaluate[D[m2[s, T, F], s]]; 
        BC1[T_?NumericQ, F_?NumericQ] := 
         NIntegrate[Cos[m1[s, T, F]], {s, 0, l}] - DD; 
        BC2[T_?NumericQ, F_?NumericQ] := -dm1[l, T, F] + dm2[l, T, F]; 
       sol2 =
          Monitor[
          FindRoot[{BC1[T, F] == 0, BC2[T, F] == 0}, {{T, 1}, {F, 1}}], {T, F}];

then define a piecewise function which is the conbination of ϕ1 and ϕ2,

fin[s_] := Piecewise[{{ϕ1[s], 0 <= s < l}, {ϕ2[s], l <= s <= L}}];

Since for different initial values used in Findroot, we can get different ϕ1 and ϕ2. However, the shape of fin[s] shown in the following figure is actually my goal.

Actually in my earlier question I have received two answers from @AlexTrounev and @bbgodfrey. However, the parameters there are set to l = L*(3/10); DD = 1; . Now in this one I especially need the results in condition of smaller DD and l. I have found there codes work not very well at small DD and l. It seems there exists a bottleneck in function FindRoot and NDSolveValue at extreme conditions, for example, in this question they are set to

l = L*(3/1000); DD = 0.01;

I’m curious about what causes such difficulty in solving these equations at this condition in function FindRoot and NDSolveValue. Are there any techniques that can help improve this situation?

For ease of reading, here I paste the answer I have received in my eariler question:

Answer from AlexTrounev:

L = 10; l = (3/10); DD = 1; mu = 1; {ϕ1''[s] - (T - F)*Sin[ϕ1[s]] == 
  0, ϕ2''[s] - T*Sin[ϕ2[s]] == 0, f1'[s] == Cos[ϕ1[s]],
  f2'[s] == Sin[ϕ2[s]], ϕ1[0] == 0, ϕ2'[L] == 
  MM, ϕ1[l] == ArcTan[1/mu], ϕ2[l] == ArcTan[1/mu], 
 f1[0] == 0, 
 f2[0] == 0}; eqs = {p1''[s]/L^2 - (T - F)*Sin[p1[s]] == 0, 
  p2''[s]/L^2 - T*Sin[p2[s]] == 0, 
  f1'[s]/L == Cos[p1[s]]}; bc = {p1[0] == 0, p1[l] == ArcTan[1/mu], 
  p2[l] == ArcTan[1/mu], f1[0] == 0, p2[1] == 0, 
  p2'[l]/L - p1'[l]/L == 0, f1[l] - DD == 0}; rule = {p1'' -> ddp1, 
  p1' -> dp1, p2'' -> ddp2, p2' -> dp2, f1' -> df1};
eqs1 = eqs /. rule; bc1 = bc /. rule;
UE[m_, t_] := EulerE[m, t];
psi[k_, n_, m_, t_] := 
  Piecewise[{{2^(k/2) UE[m, 2^k t - 2 n + 1], (n - 1)/2^(k - 1) <= t <
       n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] := 
 Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 6; M0 = 7; With[{k = k0, M = M0}, 
 nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; ycol = 
 Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 = 
 With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y;

P1 = Array[pp1, {nn}]; P2 = Array[pp2, {nn}]; F1 = Array[ff1, {nn}];
ddp1[s_] := P1 . Psi[s]; dp1[s_] := P1 . int1[s] + a0;
p1[s_] := P1 . int2[s] + a0 s + a1;
ddp2[s_] := P2 . Psi[s]; dp2[s_] := P2 . int1[s] + b0;
p2[s_] := P2 . int2[s] + b0 s + b1;
df1[s_] := F1 . Psi[s]; f1[s_] := F1 . int1[s] + c0;
eqn = Table[eqs1, {s, ycol}];
eqsAll = Join[Flatten[eqn], bc];
var = Join[{T, F, MM, a0, a1, b0, b1, c0, d0}, P1, P2, F1, F2];  
    sol2 = FindRoot[eqsAll, Table[{var[[i]], 1/10}, {i, Length[var]}]];

plot2 = Plot[Evaluate[{p1[s/L], p2[s/L]} /. sol2], {s, 0, L}, 
  PlotLegends -> {ϕ1, ϕ2}, PlotRange -> All]

fin[s_] := 
  Piecewise[{{p1[s/L] /. sol2, 0 <= s < l L}, {p2[s/L] /. sol2, 
     l L <= s <= L}}];

Plot[Evaluate[fin[s]], {s, 0, L}]

Answer from bbgodfrey:

sL = 10; sm = sL*(3/10); dd = 1; mu = 1;
sol1 = NDSolveValue[{ϕ1''[s] - g[s]*Sin[ϕ1[s]] == 0, g'[s] == 0, 
    d'[s] == Cos[ϕ1[s]], ϕ1[0] == 0, ϕ1[sm] == ArcTan[1/mu], 
    d[0] == 0, d[sm] == dd}, {ϕ1[s], ϕ1'[s], g[sm]}, {s, 0, sL}];

sol1[[3]
(* -1.13517 *)
NIntegrate[Cos[sol1[[1]]], {s, 0, sm}]
(* 1. *)
sol2 = NDSolveValue[{ϕ2''[s] - t[s]*Sin[ϕ2[s]] == 0, t'[s] == 0, 
    ϕ2'[sm] == (sol1[[2]] /. s -> sm), ϕ2[sm] == ArcTan[1/mu],
    ϕ2[sL] == 0}, {ϕ2[s], ϕ2'[s], t[sm]}, {s, 0, sL}, 
    InitialSeeding -> t[sL] == -1/2];

sol2[[3]]
(* -0.597091 *)

Answer 1

This problem has no unique solution. We can define exact solutions satisfier to ϕ1[0] ==0, ϕ2'[L] == 0 in a form

ϕ1[s_] = 2 JacobiAmplitude[k1 Sqrt[F - T] s, k1^-2];
ϕ2[s_] = Pi + 2 JacobiAmplitude[(s k2 Sqrt[T] + c), k2^-2];

where c, k1, k2 are arbitrary constants. Note that ϕ1[0] =0. Using boundary condition ϕ2'[L] == 0 we have

 ϕ2'[L]=2 k2 Sqrt[T] JacobiDN[c + k2 L Sqrt[T], 1/k2^2]
 Reduce[JacobiDN[x, 1/k2^2] == 0, x]

Solution depends on 2 integer constant

 (C[1] | C[2]) \[Element]Integers && (x == 
4 I C[2] EllipticK[1 - 1/k2^2] + 2 C[1] EllipticK[1/k2^2] - 
 InverseJacobiDN[0, 1/k^2] || 
 x == 4 I C[2] EllipticK[1 - 1/k2^2] + 2 C[1] EllipticK[1/k2^2] + 
     InverseJacobiDN[0, 1/k2^2])

Therefore

c[n1_, n2_, n3_, k_, h_] := 
  4 I n1 EllipticK[1 - 1/k^2] + 
   2 n2 EllipticK[1/k^2] + (-1)^n3 InverseJacobiDN[0, 1/k^2] - h;

where k=k2, h=k L Sqrt[T], n1,n2 are integers, and $n3=1,2$.

These multiple solutions ϕ2[s] are very sensitive to the value $l$ and DD. For example, for sL = 10; sm = sL*(1/10); dd = 1/10; mu = 1; we have from bbgodfrey code

sol1 = NDSolveValue[{\[Phi]1''[s] - g[s]*Sin[\[Phi]1[s]] == 0, 
    g'[s] == 0, 
    d'[s] == Cos[\[Phi]1[s]], \[Phi]1[0] == 0, \[Phi]1[sm] == 
     ArcTan[1/mu], d[0] == 0, d[sm] == dd}, {\[Phi]1[s], \[Phi]1'[s], 
    g[s]}, {s, 0, sL}, InitialSeeding -> {g[sm] == -1/2}]; 

Hence T-F evaluated as

sol1[[3]] /. s -> sm

(*Out[]= -13.4709*)

From our code we have same result

(T - F) /. {T -> 50.74724219864135`, F -> 64.21815429080017`}

(*Out[]= -13.4709*)

But for T we have from bbgodfrey code

sol2 = NDSolveValue[{\[Phi]2''[s] - t[s]*Sin[\[Phi]2[s]] == 0, 
    t'[s] == 0, \[Phi]2'[sm] == (sol1[[2]] /. s -> sm), \[Phi]2[sm] ==
      ArcTan[1/mu], \[Phi]2'[sL] == 0}, {\[Phi]2[s], \[Phi]2'[s], 
    t[sm]}, {s, 0, sL}, InitialSeeding -> t[sL] == 10, 
   MaxSteps -> 10^6];
sol2[[3]]
(*51.3391*)

If we define function

fin[s_] := 
  Piecewise[{{sol1[[1]], 0 <= s < sm}, {sol2[[1]], sm <= s <= sL}}];

then we can visualize sol1, sol2 as follows

{Plot[{sol1[[1]], sol2[[1]]}, {s, 0, sL}, 
  PlotLegends -> {"\[Phi]1[s]", "\[Phi]2[s]"}], 
 Plot[Evaluate[fin[s]], {s, 0, sL}]}

Note that our solution with {T -> 50.74724219864135, F -> 64.21815429080017} looks quite different