Solution of NDSolve evaluates pointwise, but won't plot (possibly due to issue with FindRoot and replacement rules)

Question

I'm trying to solve a pair of uncoupled differential equations using NDSolve, where the initial conditions are determined by FindRoot. Here's the relevant section of the code, including error message. As you can see, the result evaluates perfectly well pointwise, but fails to plot (I suspect due to some issue with hold and the order in which things are evaluated).

V0 = 10^-10;
\[Mu] = 1;
\[Omega] = 2;
amp = 0.1;
freq = 0.1;

V[\[Phi]_] := 
  V0*(Tanh[\[Phi]/\[Omega]] + amp*Sin[Tanh[\[Phi]/\[Omega]]/freq])^2;

aini = 1;
\[Phi]ini = 5;
d\[Phi]ini = -(V'[\[Phi]ini]/V[\[Phi]ini]);

sol = NDSolve[{dummy\[Phi]''[n] + 3*dummy\[Phi]'[n] - 
       1/2*(dummy\[Phi]'[n])^3 + (3 - 1/2*(dummy\[Phi]'[n])^2)*
        V'[dummy\[Phi][n]]/V[dummy\[Phi][n]] == 0, 
     dummy\[Phi][0] == \[Phi]ini, dummy\[Phi]'[0] == d\[Phi]ini}, 
    dummy\[Phi], {n, 0, 56}, MaxStepSize -> Infinity][[1, 1]];

\[Phi][n_] := dummy\[Phi][n] /. sol;
d\[Phi][n_] := Derivative[1][dummy\[Phi]][n] /. sol;
dd\[Phi][n_] := (dummy\[Phi]^\[Prime]\[Prime])[n] /. sol;

\[Epsilon]1[n_] := 1/2*d\[Phi][n]^2;
\[Epsilon]2[n_] := -(1/\[Epsilon]1[n])*d\[Phi][n]*dd\[Phi][n];
\[Epsilon]3[n_] := 
  1/(\[Epsilon]1[n]*\[Epsilon]2[n])*(dd\[Phi][n]^2 + 
      d\[Phi][n]*ddd\[Phi][n]) - \[Epsilon]2[n];

H[n_] := Sqrt[V[\[Phi][n]]/(3 - \[Epsilon]1[n])];
ascale[n_] := aini*Exp[n];

Nini[k_?NumericQ] := 
  ni /. FindRoot[
    k == 10^2*H[ni]*ascale[ni] // Rationalize[#, 0] & // 
     Evaluate, {ni, 0, 56}, WorkingPrecision -> 20];
Nfin[k_?NumericQ] := 
  nf /. FindRoot[
    k == 10^-2*H[nf]*ascale[nf] // Rationalize[#, 0] & // 
     Evaluate, {nf, 0, 56}, WorkingPrecision -> 20];

ddzoverz[n_] := 
  2 - \[Epsilon]1[n] + 3/2 \[Epsilon]2[n] + 1/4 \[Epsilon]2[n]^2 - 
   1/2 \[Epsilon]1[n]*\[Epsilon]2[n] + 
   1/2 \[Epsilon]2[n]*\[Epsilon]3[n];

eqs[k_, n_] := {(vR^\[Prime]\[Prime])[n] + (1 - \[Epsilon]1[n])*
      Derivative[1][vR][n] + (k^2/(ascale[n]^2*H[n]^2) - ddzoverz[n])*
      vR[n] == 
    0, (vI^\[Prime]\[Prime])[n] + (1 - \[Epsilon]1[n])*
      Derivative[1][vI][n] + (k^2/(ascale[n]^2*H[n]^2) - ddzoverz[n])*
      vI[n] == 0, vR[Nini[k]] == 1/Sqrt[2 k], 
   Derivative[1][vR][Nini[k]] == 0, vI[Nini[k]] == 0, 
   Derivative[1][vI][Nini[k]] == -50 Sqrt[2/k]};

vRsol[k_] := 
  NDSolve[eqs[k, n], {vR, vI}, {n, Nini[k], Nfin[k]}, 
    MaxStepSize -> Infinity][[1, 1]];
vIsol[k_] := 
  NDSolve[eqs[k, n], {vR, vI}, {n, Nini[k], Nfin[k]}, 
    MaxStepSize -> Infinity][[1, 2]];

vRk[k_?NumericQ, n_?NumericQ] := vR[n] /. vRsol[k];
vIk[k_?NumericQ, n_?NumericQ] := vI[n] /. vIsol[k];

vRk[10^11, 33]
-7.88251*10^-7
vRk[10^11, 34]
1.84531*10^-6
vRk[10^11, 35]
2.05082*10^-6

Plot[vRk[10^11, n], {n, 33, 37}]
NDSolve::dsvar: 33.00008171428572` cannot be used as a variable.

ReplaceAll::reps: {7307.25 vR[33.0001]+0.999996 (vR^\[Prime])[33.0001]+(vR^\[Prime]\[Prime])[33.0001]==0} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

ReplaceAll::reps: {7307.25 vR[33.0001]+0.999996 (vR^\[Prime])[33.0001]+(vR^\[Prime]\[Prime])[33.0001]==0.} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

NDSolve::dsvar: 33.08171436734695` cannot be used as a variable.

ReplaceAll::reps: {6206.23 vR[33.0817]+0.999996 (vR^\[Prime])[33.0817]+(vR^\[Prime]\[Prime])[33.0817]==0} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

General::stop: Further output of ReplaceAll::reps will be suppressed during this calculation.

NDSolve::dsvar: 33.163347020408175` cannot be used as a variable.

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

I'd generally define "eqs" with set rather than making it a function and using setdelay, but when I tried it this way I got the error message:

"FindRoot: The function value {-0.000544959+k} is not a list of numbers with dimensions {1} at {ni} = {0}"

I suspect my issue is the fact that everything is a function of the parameter k, so I tried a variant on this using ParametricNDSolveValue, but because k appears as part of the initial conditions, I got an error that a starting point for the variable vR couldn't be found.

Any help or workaround would be greatly appreciated!!

Answer 1

Just in case anyone runs into a similar issue, I figured out a fix! The issue was that NDSolve was setting n to a numeric value, and then trying to evaluate. So I used Block so that it would be evaluated with local n. Here is the last segment of my code after having implemented this:

vRsol[k_] := 
  Block[{n}, 
    NDSolve[eqs[k, n], {vR, vI}, {n, Nini[k], Nfin[k]}, 
     MaxStepSize -> Infinity]][[1, 1]];
vIsol[k_] := 
  Block[{n}, 
    NDSolve[eqs[k, n], {vR, vI}, {n, Nini[k], Nfin[k]}, 
     MaxStepSize -> Infinity]][[1, 2]];

vRk[k_?NumericQ, n_?NumericQ] := vR[n] /. vRsol[k];
vIk[k_?NumericQ, n_?NumericQ] := vI[n] /. vIsol[k];

Table[vRk[10^11, n], {n, 33, 37, 1}]
{-7.88037*10^-7, 1.84941*10^-6, 2.0586*10^-6, 
 6.50416*10^-7, -2.31497*10^-6}