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

Question

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

Answer 1

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.