This problem can be solved with using the Euler wavelets collocation method described here, here, and hear. First we rationalize coefficients and map solution on the unit interval, then we have
Clear["Global`*"]
testEquations = {-505/10 L0^2 x^2 a[x] phi[
x] (1 - Tanh[20 (-1 + L0 x)]) +
D[(x^2 a[
x] D[(-1 + Log[a[x]] + Tanh[20 (-1 + L0 x)]), {x}]), {x}] ==
0, 50 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) +
D[x^2 b[x] D[(-1 + Log[b[x]] +
Tanh[20 (-1 + L0 x)]), {x}], {x}] == 0,
5/100 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) +
1/1000 D[(x^2 (1 - phi[x]) D[Log[10 (1 - phi[x])], {x}]), {x}] ==
0};
Let transform system to the collocation method
testEquations /. {a'[x] -> a1[x], a''[x] -> a2[x],
b'[x] -> b1[x], b''[x] -> b2[x], phi'[x] -> phi1[x],
phi''[x] -> phi2[x]}
Out[]= {2 x a[x] (a1[x]/a[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) +
x^2 a1[x] (a1[x]/a[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) -
101/2 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) +
x^2 a[x] (-(a1[x]^2/a[x]^2) + a2[x]/a[x] -
800 L0^2 Sech[20 (-1 + L0 x)]^2 Tanh[20 (-1 + L0 x)]) == 0,
2 x b[x] (b1[x]/b[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) +
x^2 b1[x] (b1[x]/b[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) +
50 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) +
x^2 b[x] (-(b1[x]^2/b[x]^2) + b2[x]/b[x] -
800 L0^2 Sech[20 (-1 + L0 x)]^2 Tanh[20 (-1 + L0 x)]) ==
0, (-2 x phi1[x] - x^2 phi2[x])/1000 +
1/20 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) == 0}
This system we transform to the algebraic system as follows
UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] 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 = 7; M0 = 8; 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}];
xcol = 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;
var1 = Array[v1, {nn}]; var2 = Array[v2, {nn}]; var3 =
Array[v3, {nn}]; con = Array[c, {6}];
L0 = 4; a2[x_] = var1 . Psi[x]; a1[x_] = var1 . int1[x] + c[1];
a[x_] = var1 . int2[x] + c[1] x + c[2];
b2[x_] = var2 . Psi[x]; b1[x_] = var2 . int1[x] + c[3];
b[x_] = var2 . int2[x] + c[3] x + c[4];
phi2[x_] = var3 . Psi[x]; phi1[x_] = var3 . int1[x] + c[5];
phi[x_] = var3 . int2[x] + c[5] x + c[6]; var =
Join[var1, var2, var3, con];
eq = Flatten[
Table[{2 x a[x] (a1[x]/a[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) +
x^2 a1[x] (a1[x]/a[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) -
101/2 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) +
x^2 a[x] (-(a1[x]^2/a[x]^2) + a2[x]/a[x] -
800 L0^2 Sech[20 (-1 + L0 x)]^2 Tanh[20 (-1 + L0 x)]) == 0,
2 x b[x] (b1[x]/b[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) +
x^2 b1[x] (b1[x]/b[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) +
50 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) +
x^2 b[x] (-(b1[x]^2/b[x]^2) + b2[x]/b[x] -
800 L0^2 Sech[20 (-1 + L0 x)]^2 Tanh[20 (-1 + L0 x)]) ==
0, (-2 x phi1[x] - x^2 phi2[x])/1000 +
1/20 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) == 0}, {x,
xcol}]];
Finally we add boundary conditions and solve with FindRoot
bc = {a[x] == 1, b[x] == 0, phi[x] == 1} /.
x -> 1; bc1 = {a1[x] == 0, b1[x] == 0, phi1[x] == 0} /. x -> 0;
eqs = Join[eq, bc, bc1];
sol = FindRoot[eqs, Table[{var[[i]], 1/10}, {i, Length[var]}]];
Visualization
Plot[Evaluate[{a[x/L0] /. sol, b[x/L0] /. sol, phi[x/L0] /. sol}], {x,
0, L0}, PlotLegends -> {"a", "b", "phi"}]
