Plotting interpolating functions from ParametricNDSolveValue? [duplicate]

This question already has an answer here:

I'm looking at the behavior of a system of ODEs as a function of a parameter.

This is a system of phase oscillators, with the dynamics (simplified for clarity of reading, proper code farther below):

Phi'[m][t] == Freq + Sin[Phi[m + 1][t] - Phi[m][t] - p]


I want to look at the difference between two dynamic variables (for example, Phi[m+1]-Phi[m]) as a function of the parameter, p. I know how to look at the behavior of the first dynamic variable as a function of a parameter using ParametricNDSolveValue (the first Parameter Sweep problem under the Applications section of the help is similar to my approach, as is this Stack Exchange approach: Plot multiple solutions of a system of differential equations with different values of constant c by using NDsolve.)

(*system parameters*)
a = {{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {1, 0, 0, 0}};
NatFreq = {1, .98, 1.02, 1.01};
InitCond[n_] := Array[Pi (# - 1)/(2) &, n]

(*ODEs and initial conditions*)
deqns =
Table[
Phi[m]'[t] ==
NatFreq[[m]] + Sum[a[[m,j]]*Sin[Phi[j][t] - Phi[m][t] - p],{j,4}],
{m, 1, 4}]
ics = Table[Phi[m][0] == InitCond[4][[m]], {m, 1, 4}];

(*compute the solution*)
pfun1 = ParametricNDSolveValue[{deqns, ics},
Table[Phi[m][t], {m, 1, 4}], {t, 0, 20}, {p}];
F[t_, p_] := With[{val1 = First[pfun1[p][t]]},val1];


So--I am able to plot Phi[1][t] for various parameters, p, with

pVals = {0, Pi/2, Pi, 3*Pi/2};
Plot[Evaluate[Table[F[t,p], {p, pVals}]],{t, 0, 20}]


I can't figure out how to access other dynamic variables, such as Phi[2][t]. First[pfun1], best as I can tell, returns Phi[1], pfun1[[2]] is not Phi[2]. How do I plot Phi[2][t], or better yet, Phi[2][t]-Phi[1][t]?

marked as duplicate by MarcoB, Feyre, corey979, Alexey Popkov, YoungFeb 5 '17 at 17:07

• Thanks all, useful link and two helpful examples. – KBL Feb 8 '17 at 18:32

Just use this,

Plot[Evaluate[Table[pfun1[p] /. t -> x, {p, pVals}]], {x, 0, 20},PlotRange -> All,
PlotStyle -> {Black, Red, Green, Blue},Frame -> True]


For this answer credit goes to @bbgodfrey for his response to same sort of question.

a = {{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {1, 0, 0, 0}};
NatFreq = {1, .98, 1.02, 1.01};
InitCond[n_] := Array[Pi (# - 1)/(2) &, n]

deqns =
Table[
Phi[m]'[t] ==
NatFreq[[m]] + Sum[a[[m, j]]*Sin[Phi[j][t] - Phi[m][t] - p], {j, 4}],
{m, 4}];
ics = Table[Phi[m][0] == InitCond[4][[m]], {m, 4}];
vars = Table[Phi[m], {m, 4}];

pF = ParametricNDSolveValue[{deqns, ics}, vars, {t, 0, 20}, {p}];


Now to get the individual interpolation functions for a given value of p, I define the function pPhi, where the 1st argument is the Phi index and 2nd argument is the parameter value.

pPhi[n : (1 | 2 | 3 | 4), p_?NumericQ] := pF[p][[n]]


Then to plot Phi[2] for your list of p-value, I write

pVals = {0, Pi/2, Pi, 3 Pi/2};
Plot[Evaluate[Table[pPhi[2, p][t], {p, pVals}]], {t, 0, 20}]


and I get

plot