# Problem with a numeric 1/0 in NDSolveValue

The following code tries to solve three coupled differential equations with some fixed parameters

\$PrePrint = # /. {Csc[z_] :> 1/Defer@Sin[z],
Sec[z_] :> 1/Defer@Cos[z], Cot[z_] :> Defer@Cos[z]/Defer@Sin[z],
Csch[z_] :> 1/Defer@Sinh[z], Sech[z_] :> 1/Defer@Cosh[z],
Coth[z_] :> Defer@Cosh[z]/Defer@Sinh[z]} &;


These are the parameters

    ϵ := 1/100
L := 10
d := 1
ϕ := 0
χm := 5 π/12
χp := 4 π/12
αp := 1
N5 := 1
M5 := 0
N3 := 100
ΔN3 := 0
gYM := 1
δ =
1/2 Log[1/(
gYM^2 N5^2 (2 N3 - ΔN3)) (2 gYM^2 N3 N5^2 +
4 π^2 ΔN3^2 +
Sqrt[(2 gYM^2 N3 N5^2 + 4 π^2 ΔN3^2)^2 -
gYM^4 N5^4 (4 N3^2 - ΔN3^2)])];

α = -(M5/4) Cosh[δ] +
Sqrt[(π^2 N3)/gYM^2 + M5^2/16 Cosh[δ]^2];

αh = (gYM^2 α)/(4 π);


These are some functions with the parameters inside. I took these parameters as much simple I can

h1 = αp (-I α Sinh[v] -
M5/4 Log[
Tanh[(I π)/4 - (v - δ)/
2]]) + αp (I α Sinh[vb] -
M5/4 Log[Tanh[-((I π)/4) - (vb - δ)/2]]) //
FullSimplify;
h2 = αp αh (Cosh[v] + Cosh[vb]) // FullSimplify;
w = D[D[h1 h2, vb], v] // FullSimplify;
F1 = 2 h1 h2 D[h1, v] D[h1, vb] - h1^2 w // FullSimplify;
F2 = 2 h1 h2 D[h2, v] D[h2, vb] - h2^2 w;
subv = {v -> x[x2] + I y[x2], vb -> x[x2] - I y[x2]};
f42 = 2 ((F1 F2)/w^2)^(1/4) /. subv // FullSimplify // PowerExpand
(* 20 Sqrt[π] Cosh[x[x2]]^2 *)

ρ2 =
4 (F1 F2 w^2)^(1/4)/(h1 h2) /. subv /.
Sin[2 y[x2]] :> 2 Cos[y[x2]] Sin[y[x2]] // FullSimplify //
PowerExpand
(* 20 Sqrt[π] *)


This is the lagrangian and the diffrential equations involved

Lag = Factor@
FullSimplify@f42 (u'[x2]^2/u[x2]^2 + 2/
u[x2]^2) + ρ2 (x'[x2]^2 + y'[x2]^2);

equ = Factor@FullSimplify@(D[Lag, u[x2]] - D[D[Lag, u'[x2]], x2]);

eqx = Factor@FullSimplify@(D[Lag, x[x2]] - D[D[Lag, x'[x2]], x2]);

eqy = Factor@FullSimplify@(D[Lag, y[x2]] - D[D[Lag, y'[x2]], x2]);

pdes = {equ == 0, eqx == 0, eqy == 0};


These are the boundary conditions

x20 = -d Cos[ϕ];
x21 = +d Cos[ϕ];

u0 = ϵ Sqrt[1 + ((L - d Sin[ϕ])/ϵ)^2];
x0 = ArcSinh[(L - d Sin[ϕ])/ϵ];
y0 = π/2 - χm;

u1 = ϵ Sqrt[1 + ((L + d Sin[ϕ])/ϵ)^2];
x1 = ArcSinh[(L + d Sin[ϕ])/ϵ];
y1 = π/2 - χp;

Ld = L + d Sin[ϕ] // N;
L - d Sin[ϕ] > 0;

in = {u0, x0, y0} // N;
out = {u1, x1, y1} // N;

{x20, x21} // N;


Here I recollect the boundary conditions

bcs = {u[x20] == u0, x[x20] == x0, y[x20] == y0, u[x21] == u1,
x[x21] == x1, y[x21] == y1};

bcsin = {u[x20] == u0, x[x20] == x0, y[x20] == y0};

bcsd = {bcsin, u' == 0, x' == 0, y' == 0} // Flatten;


Here I try to solve the equations

sols = NDSolve[{pdes, bcs}, {u, x, y}, {x2, x20, x21}];


Here I plot the solutions

Plot[sols[][x2], {x2, x20, x21}]
Plot[sols[][x2], {x2, x20, x21}]
Plot[sols[][x2], {x2, x20, x21}]


This code produces the following errors

Power::infy: Infinite expression 1/0. encountered. >>

Power::infy: Infinite expression 1/0.^2 encountered. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

Power::infy: Infinite expression 1/0. encountered. >>

General::stop: Further output of Power::infy will be suppressed during this calculation. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

General::stop: Further output of Infinity::indet will be suppressed during this calculation. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at x2 == -1.. >>

And, when I try to plot

NDSolve::dsvar: -0.999959 cannot be used as a variable. >>

NDSolve::dsvar: -0.959143 cannot be used as a variable. >>

NDSolve::dsvar: -0.918326 cannot be used as a variable. >>

General::stop: Further output of NDSolve::dsvar will be suppressed during this calculation. >>

Notice that even the equations are 2nd order, I give only fixed points (BVP). I think the shooting method cannot be applied when I have more than one equations.

• Must be that you divide by zero... somewhere in your convoluted code. – Henrik Schumacher Mar 2 '18 at 18:58
• @HenrikSchumacher I understand that. The problem is where it could happen, because the equations involved, if you see, are quite simple. – Patrick El Pollo Mar 2 '18 at 19:02
• @corey979 Those arent variables, those are the extreme values of x2, I mean x2=[x20,x21] – Patrick El Pollo Mar 2 '18 at 19:27