Unable to find initial conditions in ParametricNDSolveValue for a system of three ODEs

Question

In my code

$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]} &;

\[Epsilon] := 1/100
L := 10
d := 1
\[Phi] := \[Pi]/4
\[Chi]m := \[Pi]/6
\[Chi]p := \[Pi]/3
\[Alpha]p := 1
N5 := 1
M5 := 1
N3 := 100
\[CapitalDelta]N3 := 4
gYM := 10
\[Delta] = 
  1/2 Log[1/(
     gYM^2 N5^2 (2 N3 - \[CapitalDelta]N3)) (2 gYM^2 N3 N5^2 + 
       4 \[Pi]^2 \[CapitalDelta]N3^2 + 
       Sqrt[(2 gYM^2 N3 N5^2 + 4 \[Pi]^2 \[CapitalDelta]N3^2)^2 - 
         gYM^4 N5^4 (4 N3^2 - \[CapitalDelta]N3^2)])];
\[Alpha] = -(M5/4) Cosh[\[Delta]] + 
   Sqrt[(\[Pi]^2 N3)/gYM^2 + M5^2/16 Cosh[\[Delta]]^2];
\[Alpha]h = (gYM^2 \[Alpha])/(4 \[Pi]);

h1 = \[Alpha]p (-I \[Alpha] Sinh[v] - 
      M5/4 Log[
        Tanh[(I \[Pi])/4 - (v - \[Delta])/
          2]]) + \[Alpha]p (I \[Alpha] Sinh[vb] - 
      M5/4 Log[Tanh[-((I \[Pi])/4) - (vb - \[Delta])/2]]);

h2 = \[Alpha]p \[Alpha]h (Cosh[v] + Cosh[vb]);
w = D[D[h1 h2, vb], v];
F1 = 2 h1 h2 D[h1, v] D[h1, vb] - h1^2 w;
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;
\[Rho]2 = 4 (F1 F2 w^2)^(1/4)/(h1 h2) /. subv;

Lag = \[Sqrt](f42/
     u[x2]^2 (f42/u[x2]^2 u'[x2]^2 + f42/
       u[x2]^2 + \[Rho]2 (x'[x2]^2 + y'[x2]^2)));
pu = D[Lag, u'[x2]];

equ = D[Lag, u[x2]] - D[D[Lag, u'[x2]], x2];

eqx = D[Lag, x[x2]] - D[D[Lag, x'[x2]], x2];

eqy = D[Lag, y[x2]] - D[D[Lag, y'[x2]], x2];

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


x20 = -d Cos[\[Phi]];
x21 = +d Cos[\[Phi]];

u0 = Sqrt[\[Epsilon]^2 + (L - d Sin[\[Phi]])^2];
x0 = ArcSinh[(L - d Sin[\[Phi]])/\[Epsilon]];
y0 = \[Pi]/2 - \[Chi]m;

u1 = Sqrt[\[Epsilon]^2 + (L + d Sin[\[Phi]])^2];
x1 = ArcSinh[(L + d Sin[\[Phi]])/\[Epsilon]];
y1 = \[Pi]/2 - \[Chi]p;

L - d Sin[\[Phi]] \[GreaterSlantEqual] 0

in = {u0, x0, y0} // N

out = {u1, x1, y1} // N

{x20, x21}

s = ParametricNDSolveValue[
  Flatten[{pdes, u[x20] == u0, u'[x20] == p, x[x20] == x0, 
    x'[x20] == q, y[x20] == y0, y'[x20] == r}], {u, x, y}, {x2, x20, 
   x21}, {p, q, r}, 
  Method -> {"EquationSimplification" -> "Residual"}]

where I try to solve a system of differential equations with free first derivatives in the initial points produces the following message

ParametricNDSolveValue::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.

when I set values for the first derivatives (parameters), for example (my goal is to fix these values in order to attain a second fixed point)

s[300, 2, -1]

I mean, the code solved "simbolically" the system, and fixing the derivatives produces the error?

You can see that I impose Method -> {"EquationSimplification" -> "Residual"} since the system is large and everything is mixed.

Looking for an answer in Documentation I got

This message is generated when NDSolve is unable to find initial conditions that satisfy the equations.

and

An error occurs because these equations do not have a solution

Notice also that this problem is related to a question I asked before in Solving a simple BVP (with an error) when only had one equation and tried to prove a simple case. The methode there worked to find the first derivative that allows to attain a second point.