Can I speed up ParametricNDSolve by converting ODE to DAE?

Question

I use ParametricNDSolve to find a solution to a very cumbersome ODE:

pf2 = ParametricNDSolveValue[{ode2$[z, k, \[Beta], 
     Rm, {q, Km}, \[CapitalLambda]] == 0, \[Phi][0] == 
    1, \[Phi]'[0] == 0}, {\[Phi]'[zR$], \[Phi][zR$]}, {z, 
   0, (1 + 10^-7) zR$}
  , {{\[Beta], 0, 1}, {\[CapitalLambda], 1, 1000}}
  , Method -> {"EquationSimplification" -> "Residual"}
  ]

Equation ode2$[z, k, \[Beta], Rm, {q, Km}, \[CapitalLambda]] is too big to show here. The only significant thing is that the parameters Rm, q, Km are given by numbers in some Do loop, and the parameters \[CapitalLambda], \[Beta] are free. The parametric function pf2[\[Beta], \[CapitalLambda]] returned by ParametricNDSolveValue is then used to construct another equation:

eqn2[\[Beta]_?NumericQ, \[CapitalLambda]_?NumericQ] = 
 pf2[\[Beta], \[CapitalLambda]] . {1, -BC1$[k, \[Beta], 
     Rm, {q, Km}, \[CapitalLambda]]}

I find the root of the equation eqn2[\[Beta], \[CapitalLambda]] ==0 given the value \[CapitalLambda]. All of this works great, but in the particular case of k==4 the calculations take about 300 times longer than with other values of k.

The reason is probably that for k==4 the ode2 equation contains 73 calls to the Appell hypergeometric function AppellF1. They occur both in the numerator and denominator of the coefficients of the equation. However, out of these 73 functions, only 7 have a different set of arguments. So I'm wondering how to make it so that I can first calculate these 7 functions and then re-use the calculated values in the coefficients of the equation. Here's how I came up with the equation modification:

"original equation ODE of 2nd order, too long to show "
ode[z, 4, p, R, \[CapitalLambda]];
"compose list of all AppellF1 functions"
n = 1;
list = {};
ode[z, 4, p, 
   R, \[CapitalLambda]] //. {AppellF1[a_, b1_, b2_, c_, x_, 
     y_] :> {AppendTo[list, AppellF1[a, b1, b2, c, x, y]];}};
"number of such functions"
list // Length
"list of unique functions and number of its occurences"
Tally[list]
list = Sort@DeleteDuplicates[list];
% // Length
"compose Rule to make substitutions in ODE"
n = 1;
rule4 = (# -> F[4, n++]) & /@ list
"compose algebraic equations for substituted AppellF1 functions"
n = 1;
setOde[z_, 4, p_, R_] = 
 Table[F[4, n], {n, Length[list]}] == list // Thread
"modifield equation ODE of 2nd order, too long to show "
ode[z, 4, p, R, \[CapitalLambda]] /. rule4;

This code gives this result:

"original equation ODE of 2nd order, too long to show "
"compose list of all AppellF1 functions"
"number of such functions"
73
"list of unique functions and number of its occurences"
{{AppellF1[1/4, -(1/2), 1/2, 5/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))], 
  25}, {AppellF1[1/4, 1/2, 1/2, 5/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))], 
  17}, {AppellF1[5/4, 1/2, 1/2, 9/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))], 
  17}, {AppellF1[5/4, -(1/2), 3/2, 9/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))], 
  9}, {AppellF1[5/4, 1/2, 3/2, 9/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))], 
  2}, {AppellF1[9/4, 1/2, 3/2, 13/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))], 
  2}, {AppellF1[9/4, -(1/2), 5/2, 13/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))], 1}}
7
"compose Rule to make substitutions in ODE"
{AppellF1[1/4, -(1/2), 1/2, 5/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))] -> F[4, 1], 
 AppellF1[1/4, 1/2, 1/2, 5/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))] -> F[4, 2], 
 AppellF1[5/4, -(1/2), 3/2, 9/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))] -> F[4, 3], 
 AppellF1[5/4, 1/2, 1/2, 9/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))] -> F[4, 4], 
 AppellF1[5/4, 1/2, 3/2, 9/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))] -> F[4, 5], 
 AppellF1[9/4, -(1/2), 5/2, 13/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))] -> F[4, 6], 
 AppellF1[9/4, 1/2, 3/2, 13/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))] -> F[4, 7]}
"compose algebraic equations for substituted AppellF1 functions"
{F[4, 1] == 
  AppellF1[1/4, -(1/2), 1/2, 5/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))], 
 F[4, 2] == 
  AppellF1[1/4, 1/2, 1/2, 5/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))], 
 F[4, 3] == 
  AppellF1[5/4, -(1/2), 3/2, 9/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))], 
 F[4, 4] == 
  AppellF1[5/4, 1/2, 1/2, 9/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))], 
 F[4, 5] == 
  AppellF1[5/4, 1/2, 3/2, 9/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))], 
 F[4, 6] == 
  AppellF1[9/4, -(1/2), 5/2, 13/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))], 
 F[4, 7] == 
  AppellF1[9/4, 1/2, 3/2, 13/
   4, -((2 p)/(-2 p + R^2)), -((2 p)/(-2 p + bv[z]^2))]}
"modifield equation ODE of 2nd order, too long to show "

Finally, we solve DAE. I'm planning on doing something like

p0Rule2 = {p0 -> (R^2 \[Beta])/(2 (-1 + R^2 + \[Beta]))};
bvRule2 = {bv[z] -> bv[z, {q, K}], 
\!\(\*SuperscriptBox[\(bv\), 
TagBox[
RowBox[{"(", "n_", ")"}],
Derivative],
MultilineFunction->None]\)[z] -> 
\!\(\*SuperscriptBox[\(bv\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", 
RowBox[{"{", 
RowBox[{"0", ",", "0"}], "}"}]}], ")"}],
Derivative],
MultilineFunction->None]\)[z, {q, K}]};
bv$[z_, {q_, K_}] = 1 + (-1 + K) Sin[(\[Pi] z)/2]^q;
ode2$[z_, k, [Beta], R, {q_, K_}, [CapitalLambda]] = 
  ode[z, k, p0, R, [CapitalLambda]] /. p0Rule2 /. bvRule2 /. 
   bv -> bv$;
setOde$[z, k, [Beta], R, {q_, K_}] = 
  setOde[z, k, p0, R] /. p0Rule2 /. bvRule2 /. bv -> bv$;
pf2 = ParametricNDSolveValue[{ode2$[z, k, [Beta], 
     Rm, {q, Km}, [CapitalLambda]] == 0, 
   setOde$[z, k, [Beta], Rm, {q, Km}], [Phi][0] == 1, [Phi]'[0] == 
    0}, {[Phi]'[zR$], [Phi][zR$]}, {z, 
   0, (1 + 10^-7) zR$}, {{[Beta], 0, 1}, {[CapitalLambda], 1, 
    1000}}, Method -> {"EquationSimplification" -> "Residual"}]

However, on this site I found many answers, the authors of which, on the contrary, advise converting the system of differential-algebraic equations DAE into the system of differential equations ODE. So I ask: What are the ways to speed up the calculations in my example?

EDIT #1: It seems that ParametricNDSolve in contrast to NDSolve cannot handle DAEs. My code suggested above does not work, and careful reading of documentation confirms that it should not work.