Return to Answer

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.

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

Plot[Evaluate[{a[x/L0] /. sol, b[x/L0] /. sol, phi[x/L0] /. sol}], {x,
   0, L0}, PlotLegends -> {"a", "b", "phi"}]

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"];

Plot[Evaluate[{a[x/L0] /. sol, b[x/L0] /. sol, phi[x/L0] /. sol}], {x,
   0, L0}, PlotLegends -> {"a", "b", "phi"}, Exclusions -> None]

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 = 5; M0 = 6; 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}]];

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