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] == 0, x'[0] == 0, y'[0] == 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[[1]][x2], {x2, x20, x21}]
Plot[sols[[2]][x2], {x2, x20, x21}]
Plot[sols[[3]][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.