# Problem of step size effectively zero

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$$.

• Are you looking solution with NDSolve only or with any other methods? – Alex Trounev Mar 2 at 17:22
• Any method that works will be fine, it is highly unlikely that this system has an analytical solution though. – Sebastián Frades Mar 2 at 22:02

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

• @SebastiánFrades You are welcome! This is more effective method for BVP then "Shooting" method implemented in NDSolve`. – Alex Trounev Mar 3 at 17:05