How to plot parametric time dependent plot with different parameters taking average of time?

Question

I am trying to plot the solution given in the code with respect to "delc". Now the problem is that it can be plot for a particular value of "t" like t=10,20,50,60,.. upto 100 but what i need is mean of all these plots. Is it possible to plot that by adjusting the y-axis so that i can get one single graph in which time also varies upto 100.

 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}, Function[t,1/2*(V11[t] + V22[t] - 2*V12[t])^(-1)], {t, 0, 
100}, delc];
 Plot[s[delc][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]

Answer 1

functions = {V11, V12, V13, V14, V21, V22, V23, V24, V31, V32, V33, 
   V34, V41, V42, V43, V44};

ClearAll[foo1]
foo1[d_, t_] := Through[s[d][t]]


delc = .1;
Plot[Evaluate@foo1[delc, t], {t, 0, 60}, Frame -> True, 
 FrameLabel -> {Style["Time", Bold, 20], Style[Subscript["S", "q"], Italic, Bold, 20]}, 
 PlotLegends -> Through[Through[functions[delc]][t]],
 ImageSize -> Large, FrameTicksStyle -> 20, PlotStyle -> Thick]

Plot[Mean@foo1[delc, t], {t, 0, 60}, Frame -> True, 
 FrameLabel -> {Style["Time", Bold, 20], Style[Subscript["S", "q"], Italic, Bold, 20]}, 
 PlotLegends -> {"Mean"}, ImageSize -> Large, FrameTicksStyle -> 20, 
 PlotStyle -> Thick]

Plot[Evaluate[foo1[delc, 50]], {delc, 0, 1}, Frame -> True, 
 FrameLabel -> {Style["Time", Bold, 20], Style[Subscript["S", "q"], Italic, Bold, 20]}, 
 PlotLegends -> LineLegend[Automatic, Through[Through[functions[ \[Delta]]][50]], 
   LegendLabel -> "t=50"], 
 ImageSize -> Large, FrameTicksStyle -> 20, PlotStyle -> Thick]

To plot 1/2*(V11[delc][t] + V22[delc][t] - 2*V12[delc][t])^(-1) you can do:

ClearAll[foo2]
foo2[d_, t_] := 1/2 /{1, 1, -2}.Through[s[d][[{1, 6, 2}]]@t];

Plot[foo2[delc,t], {t, 0, 60}, Frame -> True, 
 FrameLabel -> {Style["Time", Bold, 20], Style[Subscript["S", "q"], Italic, Bold, 20]}, 
 PlotLegends -> {Style[ToString[1/2*(V11[delc][t] + V22[delc][t] - 2*V12[delc][t])^(-1), 
     TraditionalForm], 20]}, 
 ImageSize -> Large, FrameTicksStyle -> 20, PlotStyle -> Thick]