I took the excellent answer and subsequent comment by george2079 as a challenge to obtain accurate answers to the original question in a (nearly) automated manner. It is convenient to begin by determining the value of z0
for with z''[0] == 0
, because the ODE and boundary conditions cannot be satisfied for z''[0] > 0
.
z''[0] /. Flatten@Solve[eqs[z0] /. x -> 0 /. {v[0] -> v0, z[0] -> z0, v'[0] -> 0,
z'[0] -> 0}, {v''[0], z''[0]}]
f[v0_] = (xxx /. Flatten@Solve[0 == z0 % /. z0^(17/6) -> xxx, xxx])^(6/17)
LogPlot[{f[v0], g[v0]}, {v0, -4, 4}, AxesLabel -> {v0, f}]
(* 2^(6/17) (1/(1 + Tanh[2 v0]))^(6/17) *)

Next, define the parametric solution to the ODE
s = ParametricNDSolve[{eqs[z0], v[0] == v0, z[0] == z0, v'[0] == 0, z'[0] == 0}, {v, z},
{x, 0, l/2}, {v0, z0, wp}, WorkingPrecision -> wp, MaxSteps -> Infinity];
Then, a more accurate determination of the range of z0
for which s
can be computed to at least as far as l/2
can be obtained from (here for v0 == 0
)
plt = Plot[Quiet@(z[#, z0, 15]["Domain"] /. s)[[1, 2]], {z0, 1, f[#]},
PlotRange -> All] &[0]
Flatten[{First@First[#], First@Last[#]} & /@ DeleteCases[SplitBy[
Cases[plt, Line[a__] :> a, Infinity] // Last, Last], a_ /; Length[a] == 1]]

(* {1.1521, 1.16959} *)
Thus, z0
must lie in the range just computed to solve the ODE and boundary conditions for v0 == 0
. Note that, for some value of v0
two ranges of z0
are possible. For instance, with v0 == -2
the same calculation yields

(* {2.64973, 2.8303, 5.36544, 16.8382} *)
This same computation can be performed for numerous values of v0
in the ranges {-4, 4}
in about 20 minutes on my four-processor PC.
{Range[-4, -2.4, .1], Range[-2.39, -2.31, .01], Range[-2.3, -1.6, .1],
Range[-1.59, -1.51, .01], Range[-1.5, 4, .1]} // Flatten
edges = ParallelTable[
pltt = Plot[Quiet@(z[#, z0, 15]["Domain"] /. s)[[1, 2]],
{z0, .99, f[#]}, PlotRange -> All] &[v0];
edg = {v0, Flatten[{First@First[#], First@Last[#]} & /@ DeleteCases[
SplitBy[Cases[pltt, Line[a__] :> a, Infinity] // Last, Last],
a_ /; Length[a] == 1]]},
{v0, %}]
(* {{-4., {3.84885, 283.457}}, {-3.9, {3.77187, 246.136}}, ...}*)
Plot this large table by means of
c1 = {First@#, (Last@#)[[1]]} & /@ edges;
c2 = {First@#, (Last@#)[[2]]} & /@ Cases[edges, a_ /; (Length@Last@a == 4)];
c3 = {First@#, (Last@#)[[3]]} & /@ Cases[edges, a_ /; (Length@Last@a == 4)];
c4 = {First@#, (Last@#)[[-1]]} & /@ edges;
spt = Position[edges, a_ /; (Length@a == 4)][[-1, 1]];
ListLogPlot[{c1, Join[c2, c4[[spt + 1 ;;]]], c3, c4[[;; spt]]},
PlotRange -> All, Joined -> True]

Behavior evidently is qualitatively different for v0
less than or greater than about -1
. As it turns out, l == 4
solutions for v0 > -1
are in the asymptotic regime of this nonlinear ODE. To obtain actual solutions, define
hh[v0_, bl0_, bu0_] := Module[{bl = N[bl0, 30], bu = N[bu0, 30], zt, bt},
Do[bt = (bl + bu)/2; zt = Quiet@(z[v0, bt, 30] /. s);
If[zt["Domain"][[1, 2]] < l/2 || zt[4] < 10^-4, bl = bt, bu = bt], {i, 53}]; bu]
This function takes two guesses for z0
bracketing the expected solution (i.e., one lying below the bottom curve in the plot just above, and one just below the next curve up), successively reduces the size of the bracketed region, and finally returns a quite accurate value of z0
in less than a minute. (It is based on code used to solve question 97492, which addresses a nonlinear ODE separatrix problem.) For instance, with v0 == -4
.
hh[-4, 1, 4]
Plot[Evaluate[Through[Through[{v, z}[-4, %, 30]][x]] /. s], {x, 0, 4},
PlotRange -> All]
z[-4, %%, 30][4] /. s
(* 3.79653466338964906690023326519 *)

(* 0.00010003309239505600297 *)
Thus, the boundary condition z[l/2] == 10^-4
is satisfied with a precision of better than 10^-3
. Higher precision is not difficult to obtain. Corresponding results for v0 == 4
are
hh[4, .99, f[4]]
Plot[Evaluate[Through[Through[{v, z}[4, %, 30]][x]] /. s], {x, 0, 4},
PlotRange -> All]
z[4, %%, 30][4] /. s
(* 0.999967679312725509392309112271 *)

(* 0.000100002269572818141432 *)
Although this process may seem rather lengthy, all that actually is needed to obtain an accurate numerical solution is the parametric solution s
(the second block of code above), the definition of hh
, and two reasonable initial guesses, which can be obtained as illustrated for v0 == 0
in the third block of code. Values of z0
for v0 == Range[-4, 4]
are
(* 3.79653466338964906690023326519, 3.25378157569201964527394466131,
2.64580775324334157983627578923, 1.90856596604480570876855195170,
1.15207984399450601605171562722, 1.003820075484002909681470684009,
1.000038775429121127770272366082, 0.999968958249007355677937267149
0.999967679312725509392309112271 *)
It turns out that solutions for z0
for v0
near -2
lie only between the bottom curve and the curve immediately above it in the fourth plot above.
Incidentally, I experimented with FindRoot
and FindMaximum
, which often work well with nonlinear ODE boundary-value problems, but did not obtain good results in this case.
Sample solution for l == 30
For large values of l
the entire range of v0
considered above, {-4, 4}
falls in the asymptotic regime, and obtaining solutions is somewhat more difficult. Consider l == 30
and v0 == -2
. Begin as before,
plt = Plot[Quiet@(z[#, z0, 15]["Domain"] /. ss)[[1, 2]], {z0, 1, f[#]},
PlotRange -> All] &[-2]

(Compare with third plot above.) Because the curve does not actually reach x == l/2
but instead only about x == 5
, we instead pick out the maximum point on the curve and points that bracket it. (Proceeding as we did for l = 8
would require WorkingPrecision -> 30
and a very large value of PlotPoints
, which could be prohibitively slow.)
Cases[plt, Line[a__] :> a, Infinity] // Last;
Position[%, Max[Last /@ %]][[1, 1]]
{First[%%[[% - 1]]], First[%%[[% + 1]]]}
(* 63 *)
(* {2.64485, 2.65461 *)
Now, define
dom[v0_?NumericQ, z0_?NumericQ] := Quiet@(v[v0, z0, 30]["Domain"] /. s)[[1, 2]]
gg[v0_, bl0_, bu0_] := Module[{bl = N[bl0, 30], bu = N[bu0, 30], bt, db},
Do[bt = (bl + bu)/2; db = (bu - bl)/100; If[dom[v0, bt] - dom[v0, bt - db] > 0,
bl = bt, bu = bt], {i, 68}]; bt]
Basically, it works like hh
but constructs a numerical first derivative to determine on which side of the peak in the curve that bt
lies and uses that information to reduce the size of the bracketed region on that side. Then,
rough = gg[-2, 2.645`30, 2.648`30]
z[-2, %, 30][15] /. s
(* 2.64593297120377738824924829263 *)
(* 0.588648194538264076985120037 *)
gives a rough solution for z0
. Refine it with
hh[-2, rough - 10^-15, rough]
where the Do
limit in hh
has been increased to 70
.
Plot[Evaluate[Through[Through[{v, z}[-2, %, 30]][x]] /. s], {x, 0, 15},
PlotRange -> All]
z[-2, %%, 30][15] /. s
rough - %%%
(* 2.64593297120377738824924570875 *)

(* 0.0001000074744809224181 *)
(* -2.58388*10^-24 *)
That the accurate solution differs from the rough
solution by only about one part in 10^24
is indicative of the sensitivity of such separatrix-related computations.