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]}],Jacobian -> "FiniteDifference"];
Visualization
Plot[Evaluate[{a[x/L0] /. sol, b[x/L0] /. sol, phi[x/L0] /. sol}], {x,
0, L0}, PlotLegends -> {"a", "b", "phi"}, Exclusions -> None]
We also can solve this problem with linear FEM, using iterative false transient algorithm discussed here, here, and here. We use transformed system of equations with adding iteration parameter dt=1/5 and linearized as follows
Clear["Global`*"]
Needs["NDSolve`FEM`"]
xmesh = ToElementMesh[ImplicitRegion[0 <= x <= 1, {x}],
MaxCellMeasure -> 2 10^-3]
eqn = {-(a[x] - A[i - 1][x])/dt -
808 a[x] P[i - 1][x] (1 - Tanh[20 (-1 + 4 x)]) +
2/x (80 Sech[20 (-1 + 4 x)]^2 a[x] +
Derivative[1][a][x]) + (80 Sech[20 (-1 + 4 x)]^2 Derivative[
1][a][x] +
A1[i - 1][x] a'[x]/A[i - 1][x]) + (-12800 Sech[
20 (-1 + 4 x)]^2 Tanh[20 (-1 + 4 x)] a[x] -
A1[i - 1][x] a'[x]/A[i - 1][x] + (a^\[Prime]\[Prime])[x]) ==
NeumannValue[0, x == 0], -(b[x] - B[i - 1][x])/dt +
800 P[i - 1][x] a[x] (1 - Tanh[20 (-1 + 4 x)]) +
2/x (80 Sech[20 (-1 + 4 x)]^2 b[x] +
Derivative[1][b][x]) + (80 Sech[20 (-1 + 4 x)]^2 Derivative[1][
b][x] + B1[i - 1][x] b'[x]/B[i - 1][x]) + (-12800 Sech[
20 (-1 + 4 x)]^2 Tanh[20 (-1 + 4 x)] b[x] -
B1[i - 1][x] b'[x]/B[i - 1][x] + (b^\[Prime]\[Prime])[x]) ==
NeumannValue[0, x == 0], (phi[x] - P[i - 1][x])/dt +
4/5 P[i - 1][x] a[
x] (1 - Tanh[20 (-1 + 4 x)]) + (-2/x Derivative[1][phi][x] - (
phi^\[Prime]\[Prime])[x])/1000 == NeumannValue[0, x == 0]}; bc =
DirichletCondition[{a[x] == 1, b[x] == 0, phi[x] == 1}, x == 1];
Initial guess is very important in this case, but we use the simple one
A[0][x_] := 1; B[0][x_] = 1/100; P[0][x_] = 1; A1[0][x_] := 0;
B1[0][x_] = 0; P1[0][x_] := 0;
L0 = 4; dt =
1/5; nn = 231; Do[{A[i], B[i], P[i], A1[i], B1[i], P1[i]} =
NDSolveValue[{eqn, bc}, {a, b, phi, a', b', phi'},
Element[{x}, xmesh]]; , {i, 1, nn}]
Visualization
Table[Plot[Evaluate[{A[i][x/L0], B[i][x/L0], P[i][x/L0]}], {x, 0, L0},
PlotLegends -> {"a", "b", "phi"}, PlotRange -> All,
PlotLabel -> i], {i, 228, 231}]
Also FEM solution looks similar to above computed with wavelets and computed with pdetoae at xzczd post, the algorithm not stable and solution diverges at nn>231.
