Plotting steady state behavior as a function of a parameter

Question

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]

Thanks to helpful answers here (Plotting interpolating functions from ParametricNDSolveValue?), I understand how to plot the difference between two dynamic variables (for example, dPhi = Phi[m+1]-Phi[m]) for different values of the parameter, p.

I'm still struggling to plot dPhi as a function of parameter, p. This system has a transient before reaching a steady state, and I'd like to plot this steady state phase difference versus the parameter p, either as a Plot or a ListPlot.

(*system parameters*)
a = RotateLeft[IdentityMatrix[4], {1, 0}];
NatFreq = {1, .98, 1.02, 1.01};
InitCond[n_] := Array[Pi (# - 1)/(4) &, 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}];
vars = Table[Phi[m], {m, 4}];
pVals = {0, Pi/2, Pi, 3 Pi/2};

(*compute the solution*)
pfun = ParametricNDSolveValue[{deqns, ics}, vars, {t, 0, 100}, {p}];
pPhi[n : (1 | 2 | 3 | 4), p_] := pfun[p][[n]]

With this, I can make a lovely plot of the transient of the phase differences for several values of p:

d41 = Mod[
Evaluate[
Table[pPhi[4,p][t], {p, pVals}] - 
Table[pPhi[1,p][t], {p, pVals}]
], 
2*Pi];

Plot[d41, {t, 0, 100},
PlotRange -> {All, {0, 2*Pi}},
PlotStyle -> {Black, Red, Green, Blue},
AxesLabel -> {time, Phi4 - Phi1},
PlotLegends -> pVals]

What I haven't been able to figure out is how to plot the final values of these transients against the parameter p rather than against time. This was my best guess:

T=50;
Plot[Evaluate[pPhi[2,p][t] /. t -> T-pPhi[1,p][t] /. t -> T], {p, 0, 2*Pi},
PlotRange -> {{0, 2*Pi}, {0, 2*Pi}]

Everything I try produces many errors. This seems like it may be very simple, but I haven't been able to work my way through it. (Please also note that it must be NDSolve related and not DSolve; I use Sine in this is example, but I sometimes use a longer Fourier series.) Thanks very much.

Answer 1

I often find it is easier to work with the solutions if I solve for expressions rather than just the functions. In the code below, I solved for functions of t, and then used replacement to obtain values for t = 100.

(* Solve for functions of t *)
vars = Table[Phi[m][t], {m, 4}];

(*compute the solution*)
pfun = ParametricNDSolveValue[{deqns, ics}, vars, {t, 0, 100}, {p}];

(* Plot for p= Pi *)
Plot[pfun[Pi][[4]] - pfun[Pi][[1]], {t, 0, 100}, PlotRange -> All]

(* A plot for t=100 vs p *)
Plot[(pfun[p][[4]] - pfun[p][[1]]) /. t -> 100, {p, 0, 3 Pi/2}]

Answer 2

As suggested by @march, you can do it like this,

T=50;
Plot[Mod[Evaluate[pPhi[4, p][t] /. t -> T] - Evaluate[pPhi[1, p][t] /. t -> T], 2 π], 
{p, 0, 2 π},PlotRange -> All]