# How to plot solution of first order time dependent differential equations with other parameter not with time

I am trying to plot the solution given in the code with respect to "delc", not w.r.t time "t". I don't know what to do for getting the plot between solution and value "delc" which vary 0 to 2. The initial condition are also included in the text. The only thing that i want is, plot between solution and "delc" for any arbitrary value of "t". If anyone can solve this is most welcome.

w1 = 1;
gma1 = 0.005;
n1 = 1;
gma2 = 0.005;
G1 = 0.005;
k1 = .1;
k2 = 0.1;
a1 = 0.07;
a2 = 0.58;
k0 = 0.1;
Q1 = 1.268;
del0 = 1;
N1 = 1;
ome = 1;
M1 = del0*(1 - Cos[ome*t]);
s = ParametricNDSolveValue[{V11'[t] - V21[t]*w1 - V12[t]*w1 == 0,
V12'[t] - V22[t]*w1 + w1*V11[t] + gma1*V12[t] -
Sqrt[2]*G1*a1*V13[t] - Sqrt[2]*G1*a2*V14[t] == 0,
V13'[t] - V23[t]*w1 + k1*V13[t] +
Sqrt[2]*G1*a2*V11[t] - (-G1*Q1 + delc)*V14[t] == 0,
V14'[t] - V24[t]*w1 + k1*V14[t] -
Sqrt[2]*G1*a1*V11[t] - (G1*Q1 - delc)*V13[t] == 0,
V21'[t] + V11[t]*w1 + gma1*V21[t] - Sqrt[2]*G1*a1*V31[t] -
Sqrt[2]*G1*a2*V41[t] - w1*V22[t] == 0,
V22'[t] + V12[t]*w1 + gma1*V22[t] - Sqrt[2]*G1*a1*V32[t] -
Sqrt[2]*G1*a2*V42[t] + w1*V21[t] + gma1*V22[t] -
Sqrt[2]*G1*a1*V23[t] - Sqrt[2]*G1*a2*V24[t] - gma2*(2*n1 + 1) ==
0, V23'[t] + w1*V13[t] + gma1*V23[t] - Sqrt[2]*G1*a1*V33[t] -
Sqrt[2]*G1*a2*V43[t] + k1*V23[t] +
Sqrt[2]*G1*a2*V21[t] - (-G1*Q1 + delc)*V24[t] == 0,
V24'[t] + V14[t]*w1 + gma1*V24[t] - Sqrt[2]*G1*a1*V34[t] -
Sqrt[2]*G1*a2*V44[t] + k1*V24[t] -
Sqrt[2]*G1*a1*V21[t] - (G1*Q1 - delc)*V23[t] == 0,
V31'[t] + k1*V31[t] +
Sqrt[2]*G1*a2*V11[t] - (-G1*Q1 + delc)*V41[t] - w1*V32[t] == 0,
V32'[t] + k1*V32[t] +
Sqrt[2]*G1*a2*V12[t] - (-G1*Q1 + delc)*V42[t] + w1*V31[t] -
Sqrt[2]*G1*a2*V34[t] - Sqrt[2]*G1*a1*V33[t] + gma1*V32[t] == 0,
V33'[t] + k1*V33[t] +
Sqrt[2]*G1*a2*V13[t] - (-G1*Q1 + delc)*V43[t] + k1*V33[t] +
Sqrt[2]*G1*a2*V31[t] - (-G1*Q1 + delc)*V34[t] - k0 == 0,
V34'[t] + k1*V34[t] +
Sqrt[2]*G1*a2*V14[t] - (-G1*Q1 + delc)*V44[t] + k1*V34[t] -
Sqrt[2]*G1*a1*V31[t] - (G1*Q1 - delc)*V33[t] == 0,
V41'[t] + k1*V41[t] -
Sqrt[2]*G1*a1*V11[t] - (G1*Q1 - delc)*V31[t] - w1*V42[t] == 0,
V42'[t] + k1*V42[t] +
Sqrt[2]*G1*a1*V12[t] - (G1*Q1 - delc)*V32[t] + w1*V41[t] -
Sqrt[2]*G1*a2*V44[t] - Sqrt[2]*G1*a1*V43[t] + gma1*V42[t] == 0,
V43'[t] + k1*V43[t] -
Sqrt[2]*G1*a1*V13[t] - (G1*Q1 - delc)*V33[t] + k1*V43[t] +
Sqrt[2]*G1*a2*V41[t] - (-G1*Q1 + delc)*V44[t] == 0,
V44'[t] + k1*V44[t] -
Sqrt[2]*G1*a1*V14[t] - (G1*Q1 - delc)*V34[t] + k1*V44[t] -
Sqrt[2]*G1*a1*V41[t] - (G1*Q1 - delc)*V43[t] - k0 == 0,
V11[0] == 1, V12[0] == 1, V13[0] == 0, V14[0] == 0, V21[0] == 0,
V22[0] == 1, V23[0] == 0, V24[0] == 0, V31[0] == 0, V32[0] == 0,
V33[0] == 0, V34[0] == 0, V41[0] == 0, V42[0] == 0, V43[0] == 0,
V44[0] == 0}, {V11, V12, V13, V14, V21, V22, V23, V24, V31, V32,
V33, V34, V41, V42, V43, V44},{t, 0, 100},g0];
P2 = Plot[{Evaluate[1/2*(V11[t] + V22[t] - 2*V12[t])^(-1) /. s]}, {t,
0, 60}, PlotRange -> {0, 1}, Frame -> True,
FrameLabel -> {Style["Time", Bold, 20],
Style[" \!$$\*SubscriptBox[\(S$$, $$q$$]\)", Bold, 20]},
FrameTicksStyle -> Directive[FontSize -> 20],
PlotStyle -> {Thickness[0.0005], Thickness[0.008]}]
• Can't find $g_0$ in your code for the DEs. Jan 21, 2020 at 11:34
• g0 is the parameter inside ParametricNDSolveValue . Jan 21, 2020 at 11:35
• It's unclear what you want to plot "about g0". P2=... is a time plot (which doesn't work) Jan 21, 2020 at 11:37
• @LouisB , sorry its not g0 , it is delc. Now please again see the text.
– vini
Jan 21, 2020 at 11:46

To plot the solution for a range of delc, make these changes to your code:

(1) Remove delc = 1; Leave delc undefined. It will be your parameter.

(2) Keep the ParametricDSolveValue command and its first argument, but change the other arguments to get

s = ParametricNDSolveValue[  ...   ,
{1/2*(V11[t] + V22[t] - 2*V12[t])^(-1)}, {t, 0, 100}, delc];

Notice that the second argument is just the expression that we want to plot. This is not the only way to do it, but it has its advantages. Evaluate the code and plot the solution over the range of delc at t = 60 like this

Plot[s[delc] /. t -> 60, {delc, 0, 2},
Frame -> True, PlotRange -> {All, {0.545, .5475}},
FrameLabel -> {Style["\[CapitalDelta]c", Bold, 20],
Style["\!$$\*SubscriptBox[\(S$$, $$q$$]\)", Bold, 20],
Style["t = 60", Bold, 20]},
FrameTicksStyle -> Directive[FontSize -> 20],
GridLines->Automatic]

We can also plot the solution versus time for delc = 1.25 like this

Plot[s[1.25], {t, 0, 100},
PlotRange -> {All, {0, 15}}, Frame -> True,
FrameLabel -> {Style["t", Bold, 20],
Style["\!$$\*SubscriptBox[\(S$$, $$q$$]\)", Bold, 20],
Style["\[CapitalDelta]c = 1.25", Bold, 20]},
FrameTicksStyle -> Directive[FontSize -> 20],
GridLines -> Automatic]

Edit to Plot an Average

One way to average $$S_q(\Delta c)$$ for $$\Delta c$$ range from 0 to 2 is to simply evaluate s[delc] at various values of $$\Delta c$$ and find their mean. Here is the (complete) code that accomplishes that:

ClearAll["Global`*"]
w1 = 1;
gma1 = 0.005;
n1 = 1;
gma2 = 0.005;
G1 = 0.005;

k1 = .1;
k2 = 0.1;
a1 = 0.07;
a2 = 0.58;
k0 = 0.1;
Q1 = 1.268;
del0 = 1;
N1 = 1;
ome = 1;
M1 = del0*(1 - Cos[ome*t]);

s = ParametricNDSolveValue[{
V11'[t] - V21[t]*w1 - V12[t]*w1 == 0,
V12'[t] - V22[t]*w1 + w1*V11[t] + gma1*V12[t] -
Sqrt[2]*G1*a1*V13[t] - Sqrt[2]*G1*a2*V14[t] == 0,
V13'[t] - V23[t]*w1 + k1*V13[t] +
Sqrt[2]*G1*a2*V11[t] - (-G1*Q1 + delc)*V14[t] == 0,
V14'[t] - V24[t]*w1 + k1*V14[t] -
Sqrt[2]*G1*a1*V11[t] - (G1*Q1 - delc)*V13[t] == 0,
V21'[t] + V11[t]*w1 + gma1*V21[t] -
Sqrt[2]*G1*a1*V31[t] - Sqrt[2]*G1*a2*V41[t] - w1*V22[t] == 0,
V22'[t] + V12[t]*w1 + gma1*V22[t] -
Sqrt[2]*G1*a1*V32[t] - Sqrt[2]*G1*a2*V42[t] +
w1*V21[t] + gma1*V22[t] - Sqrt[2]*G1*a1*V23[t] -
Sqrt[2]*G1*a2*V24[t] - gma2*(2*n1 + 1) == 0,
V23'[t] + w1*V13[t] + gma1*V23[t] -
Sqrt[2]*G1*a1*V33[t] - Sqrt[2]*G1*a2*V43[t] + k1*V23[t] +
Sqrt[2]*G1*a2*V21[t] - (-G1*Q1 + delc)*V24[t] == 0,
V24'[t] + V14[t]*w1 + gma1*V24[t] - Sqrt[2]*G1*a1*V34[t] -
Sqrt[2]*G1*a2*V44[t] + k1*V24[t] -
Sqrt[2]*G1*a1*V21[t] - (G1*Q1 - delc)*V23[t] == 0,
V31'[t] + k1*V31[t] + Sqrt[2]*G1*a2*V11[t] -
(-G1*Q1 + delc)*V41[t] - w1*V32[t] == 0,
V32'[t] + k1*V32[t] + Sqrt[2]*G1*a2*V12[t] -
(-G1*Q1 + delc)*V42[t] + w1*V31[t] -
Sqrt[2]*G1*a2*V34[t] - Sqrt[2]*G1*a1*V33[t] + gma1*V32[t] == 0,
V33'[t] + k1*V33[t] + Sqrt[2]*G1*a2*V13[t] -
(-G1*Q1 + delc)*V43[t] + k1*V33[t] +
Sqrt[2]*G1*a2*V31[t] - (-G1*Q1 + delc)*V34[t] - k0 == 0,
V34'[t] + k1*V34[t] + Sqrt[2]*G1*a2*V14[t] -
(-G1*Q1 + delc)*V44[t] + k1*V34[t] -
Sqrt[2]*G1*a1*V31[t] - (G1*Q1 - delc)*V33[t] == 0,
V41'[t] + k1*V41[t] - Sqrt[2]*G1*a1*V11[t] -
(G1*Q1 - delc)*V31[t] - w1*V42[t] == 0,
V42'[t] + k1*V42[t] + Sqrt[2]*G1*a1*V12[t] -
(G1*Q1 - delc)*V32[t] + w1*V41[t] -
Sqrt[2]*G1*a2*V44[t] - Sqrt[2]*G1*a1*V43[t] + gma1*V42[t] == 0,
V43'[t] + k1*V43[t] - Sqrt[2]*G1*a1*V13[t] -
(G1*Q1 - delc)*V33[t] + k1*V43[t] +
Sqrt[2]*G1*a2*V41[t] - (-G1*Q1 + delc)*V44[t] == 0,
V44'[t] + k1*V44[t] - Sqrt[2]*G1*a1*V14[t] -
(G1*Q1 - delc)*V34[t] + k1*V44[t] -
Sqrt[2]*G1*a1*V41[t] - (G1*Q1 - delc)*V43[t] - k0 == 0,
V11[0] == 1, V12[0] == 1, V13[0] == 0, V14[0] == 0,
V21[0] == 0, V22[0] == 1, V23[0] == 0, V24[0] == 0,
V31[0] == 0, V32[0] == 0, V33[0] == 0, V34[0] == 0,
V41[0] == 0, V42[0] == 0, V43[0] == 0, V44[0] == 0},
1/2*(V11[t] + V22[t] - 2*V12[t])^(-1), {t, 0, 100}, delc];

With[{avg = Mean@Table[s[dc], {dc, Subdivide[0, 2, 10]}]},
Plot[avg, {t, 0, 100}, PlotRange -> {0, 15}]]

Note that Subdivide[0, 2, 10] gives 11 evenly spaced sample points. You can generate that same list using Range[0, 2, 2/10] or any other method you choose. Eleven is enough points for the plot. It turns out s[delc] gives practically the same function of $$t$$ for any value of delc in the interval from 0 to 2.

• Please clarify 2nd point you mentioned above because its still not giving any plot.
– vini
Jan 22, 2020 at 6:44
• Here's what I did: Start with a fresh kernel. Copy / paste the code from the question into a new notebook. Remove the line containing delc = 1;. Remove all the code beginning at {V11, V12, V13 .. to the end. Copy / paste the second line of code from item (2) in answer. Copy / paste the plot command from the answer. Runs fine for me. If you can't get it to work, maybe you should edit your question to show exactly what code you are using and what results or error message you are getting. Jan 22, 2020 at 6:59
• Thanks very much. its working now.
– vini
Jan 22, 2020 at 7:07
• This plot that you have plotted, is for a particular value of time.i.e t=60, can you tell me how can i plot for t=30,10,20,70,.... and after that average of all plots?
– vini
Jan 30, 2020 at 14:34
• Plot[s[delc] /. t -> {10, 20, 30, 70}, {delc, 0, 2}] gives a quick plot at the 4 times. The "curves" look like straight lines because of the scale on the vertical axis. I don't think you really want to average these curves. Please consider asking a new question and tell us what you are trying to achieve, show or study, so we can provide answers that better hit the target. You could average those curves, however, with Plot[Mean[s[delc] /. t -> {10, 20, 30, 70}], {delc, 0, 2}] Jan 30, 2020 at 22:21