Problem of step size effectively zero

Question

I've been trying to solve the next system of differential equations which is very similar to this one in which I also sought help

Step size is effectively zero

$$F^2-G^2+HF'-F''+1=0$$

$$2GF+HG'-G''=0$$

$$2F+H'=0$$

$$HH'+P'-H''=0$$

With boundary conditions

$$F(0)=G(0)=H(0)=P(0)=0$$

$$F(∞)=0, G(∞)=1$$

Again, I turned them into a system of first order ordinary differential equations, $F=x, F'=y, G=z, G'=s, H'=p$ and $P=u$. So I tried the same code that run with the original system, with a few corrections, but it doesn't seem to work with this system, probably because it is unstable too. At t = 3.45 there is a step size effectively zero problem. I also tried with ParametricNDSolve, but I really don't know how it works. In the original system I had an idea of the behavior of the solution unlike this case in which I don't, so I assumed P(14) = 0, in an attempt to copy the idea of getting a similar solution and then improve it.

odes = {x'[t] == y[t], y'[t] == x[t]^2 - z[t]^2 + s[t] y[t] + 1, 
   s'[t] == 2 z[t] x[t] + p[t] s[t], z'[t] == s[t], p'[t] == -2 x[t], 
   u'[t] == 2 x[t] p[t] - 2 y[t]};
bcs = {x[0] == 0, z[0] == 0, p[0] == 0, x[14] == 0, z[14] == 1, 
   u[0] == 0};
vars = {x, y, z, s, p, u};
sol32 = NDSolve[{odes, bcs}, vars, {t, 0, 14}, 
   Method -> {"Shooting", 
     "StartingInitialConditions" -> {s[14] == 0, y[14] == 0, 
       p[14] == 0, x[14] == 0, z[14] == 1, u[14] == 0}}, 
   PrecisionGoal -> 6, AccuracyGoal -> 6, WorkingPrecision -> 32];

solMP = NDSolve[{odes, bcs}, vars, {t, 0, 14}, 
  Method -> {"Shooting", 
    "StartingInitialConditions" -> {Through[vars[14]] == 
       Through[vars[14] /. First[sol32]]}}, PrecisionGoal -> 10, 
  AccuracyGoal -> 10]

ListLinePlot[{x, z, p} /. First[sol], PlotLegends -> {F, G, H}]

It is the exact same one in the original post but with the equation $u'= 2 x p - 2y$, and $u(14) = 0$.

Answer 1

This problem can be solved with collocation method using Bernoulli wavelets. First we map interval to the unit interval and define wavelets, functions and derivatives for vars = {x, y, z, s, p, u} on the unit domain as follows

Clear["Global`*"]

n = 4;
M = Sum[1, {j, 0, n, 1}, {i, 0, 2^j - 1, 1}] + 1;
dx = 1/M; A = 0; xl = Table[A + l*dx, {l, 0, M}]; tcol = 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, M + 1}]; 
psi1[x_] := Piecewise[{{BernoulliB[2, x], 0 <= x < 1}, {0, True}}];
psi2[x_] := Piecewise[{{BernoulliB[1, x], 0 <= x < 1}, {0, True}}];
psi1jk[x_, j_, k_] := psi1[j*x - k];
psi2jk[x_, j_, k_] := psi2[j*x - k];
psijk = Compile[{{x, _Real}, {j, _Integer}, {k, _Integer}}, (psi1jk[x,
        j, k] + psi2jk[x, j, k])/2];
(*intjk=Integrate[psijk[x,j,k],x,Assumptions\[Rule]{j>0,k>0}]*)
psijk1 = Compile[{{x, _Real}, {j, _Integer}, {k, _Integer}}, 
   Piecewise[{{(-k + k^2)/(2*j), 
      j > 0 && k == 0 && 
       1/j - x < 0}, {(1/6)*(-x + 3*k^2*x - 3*j*k*x^2 + j^2*x^3), 
      j > 0 && k == 0 && x > 0 && 
       1/j - x >= 
        0}, {(k - k^3 - j*x + 3*j*k^2*x - 3*j^2*k*x^2 + j^3*x^3)/(6*
         j), j > 0 && k > 0 && k/j - x < 0 && 1/j + k/j - x >= 0}}, 
    0]];
Psi[x_] := 
  Join[{1}, 
   Flatten[Table[psijk[x, 2^j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]];
int1[x_] := 
  Join[{x}, 
   Flatten[Table[
     psijk1[x, 2^j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]];
var1 = Join[{a0}, 
  Flatten[Table[a[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]]; var2 = 
 Join[{b0}, 
  Flatten[Table[b[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]]; var3 = 
 Join[{c0}, 
  Flatten[Table[c[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]]; var4 = 
 Join[{d0}, 
  Flatten[Table[d[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]]; var5 = 
 Join[{e0}, 
  Flatten[Table[e[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]]; var6 = 
 Join[{f0}, Flatten[Table[f[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]];
z1[t_] := var3 . Psi[t]; z[t_] := var3 . int1[t] + c1;
y1[t_] := var1 . Psi[t]; y[t_] := var1 . int1[t] + a1; 
x1[t_] := var2 . Psi[t]; x[t_] := var2 . int1[t] + b1;
s1[t_] := var4 . Psi[t]; s[t_] := var4 . int1[t] + d1;
p1[t_] := var5 . Psi[t]; p[t_] := var5 . int1[t] + e1; 
u1[t_] := var6 . Psi[t]; u[t_] := var6 . int1[t] + f1;

Second, we define boundary conditions and equations to be optimize

L = 14;
    bcs = {x[0] == 0, z[0] == 0, p[0] == 0, x[1] == 0, z[1] == 1, 
       u[0] == 0};
  

  eq = Flatten[Table[{-x1[xcol[[j]]]/L + y[xcol[[j]]], 1 + x[xcol[[j]]]^2 + s[xcol[[j]]]*y[xcol[[j]]] - y1[xcol[[j]]]/L - z[xcol[[j]]]^2, 
      p[xcol[[j]]]*s[xcol[[j]]] - s1[xcol[[j]]]/L + 2*x[xcol[[j]]]*z[xcol[[j]]], s[xcol[[j]]] - z1[xcol[[j]]]/L, 
      -p1[xcol[[j]]]/L - 2*x[xcol[[j]]], -u1[xcol[[j]]]/L + 2*p[xcol[[j]]]*x[xcol[[j]]] - 2*y[xcol[[j]]]}, {j, M}]]; 
varM = Join[{a1, b1, c1, d1, e1, f1}, var1, var2, var3, var4, var5, var6]; 

Finally we use NMinimize to solve this problem

sol = NMinimize[{Norm[eq], bcs}, varM]

With this bcs we have norm of equations of 3.83121*10^-6 and it is not bad. For visualization we use Plot so Compile complains, it is why we use Plot with Quite

Plot[Evaluate[{x[t/L], z[t/L], p[t/L]} /. sol[[2]]], {t, 0, 14}, 
  PlotLegends -> {"F", "G", "H"}, Frame -> True] // Quiet

We can improve optimal numerical solution using FindRoot as follows

sol1 = FindRoot[Join[Table[eq[[i]] == 0, {i, Length[eq]}], bcs], 
  Table[{varM[[i]], varM[[i]] /. sol[[2]]}, {i, Length[varM]}], 
  MaxIterations -> 1000] 

But it produces same plot

Plot[Evaluate[{x[t/L], z[t/L], p[t/L]} /. sol1], {t, 0, 14}, 
  PlotLegends -> {"F", "G", "H"}, Frame -> True] // Quiet 

For some implementation of the colocation method with Bernoulli wavelets see my post on this forum

High precision numerical solution of the nonlinear Volterra integral equation

There are also answer here with application Bernoulli wavelets to optimization problem. Some papers about this method

A new numerical method for fractional order Volterra integro-differential equations

The Bernoulli wavelets operational matrix of integration and its applications for the solution of linear and nonlinear problems in calculus of variations