# Plotting steady state behavior as a function of a parameter

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.

• I'm having trouble reproducing your graph. I copied and pasted your code, ran it, and got this. Can you check your code? – march Feb 8 '17 at 23:55
• Hi, the edit made to my code by David is incorrect (I am very new to Stack Exchange, so not sure what I should do to fix this)--my matrix is not the identity matrix, it is the identity matrix shifted left by one column. In physical terms, this means that element 1 receives input from element 2, and so forth. the original a matrix is: a = {{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {1, 0, 0, 0}}; – KBL Feb 9 '17 at 0:10
• I have "rolled back" the edit, I assume this is the proper change to make! – KBL Feb 9 '17 at 0:12
• In the interests of clarity, I changed my matrix a so it is expressed in terms of the Identity Matrix, but is still correct. – KBL Feb 9 '17 at 0:17
• I think with your code, the straight-forward thing works: Plot[Mod[pPhi[4, p][100] - pPhi[1, p][100], 2 \[Pi]], {p, 0, 2 \[Pi]}]. – march Feb 9 '17 at 6:31

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


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]