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

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!!

• "ascale"and "H" and "epsiloni" are not defined. May 6, 2022 at 8:41
• Sorry--I've edited the question to include all definitions now. May 6, 2022 at 14:13
• Derivatives are not written; "vR^[Prime][Prime][n]", but "vR''[n]". Or the long form: "Derivative[2][vR][n]" May 6, 2022 at 15:22
• I'm not sure what's up with the formatting here, but in the actual notebook they are written as "vR''[n]" May 6, 2022 at 19:48
• \[Prime] has character code 8242, however the character used in "vR''[n]" is a single quote with character code 39 May 6, 2022 at 20:07

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}
`