Changing Plot Lines into Different Colors When Plotting Table of NDSolve Solutions

Question

I am using NDSolve to solve two differential equations and then putting the solutions into a table and varying one of the parameters to see a few different cases.

I am able to get everything to work and output a plot with different lines for each case, but I cannot figure out how to get these lines to be different colors and how to label them individually or at least add a legend. It seems to be viewing them as a single solution when I plot. If I try to label them using PlotLabel, it just puts the first one on the first solution and nothing on the other lines.

The plots in question are the three at the very bottom

ClearAll["Global`*"]

Q =  100 ;(*Heat into the shroud in Watts. Based on roughly 1350 W/m^2 from the solar simulator on one face of the shroud*)
QAmb = 0  ;(*Heat loss to ambinet. Zero for now*)
A =  2*1.935*10^-5 ;(*Area of orifice m^2  Based on 1/4 inch pipe with 0.028 inch wall thickness (two outlets)*)
h =  199000 ;(*Heat of vaporization of LN2 J/kg*)
Cp =  1039 ;(*Specific heat of nitrogen J/(kg K) *)
R =  296.8 ;(*Gas constant for nitrogen J/(kg K)*)
\[Gamma] = 1.40  ;(*Specific heat ratio*)
V = 0.001 ; (*Enclosed volume m^3*)
Pe = 101000 ; (*External pressure in pa*)
\[Rho]o = 4 ;(*Approx density of nitrogen at 80K in kg/m^3. This was the lowest temp data I could find*)
mi = \[Rho]o*V; (*inital mass in the chamber*)
tf = 120 ;(*Final time in seconds*)

P = m[t]*R*T[t]/V ;(*Pressure term*)
mDotEvap = Q/h ; (*rate of evap*)
mDotOut = (P*A/Sqrt[T[t]])*Sqrt[(2\[Gamma]/(R(\[Gamma]-1)))*((Pe/P)^(2/\[Gamma])-(Pe/P)^((\[Gamma]+1)/\[Gamma]))];  (*mass flow out of the orifice*)
mDotOutChoked =  (P*A/Sqrt[T[t]])*Sqrt[\[Gamma]/R]*(2/(\[Gamma]+1))^((\[Gamma]+1)/(2(\[Gamma]-1))); (*mass flow out of the orifice if choked*)

TEq = T'[t] == 1/(m[t]*Cp)(mDotEvap*h-mDotOut*Cp*T[t]-QAmb) ; (*Diff Eq for Temperature in the cavity*)
mEq = m'[t] == mDotEvap - mDotOut ; (*Conservation of mass*)

icT = T[0] == 77 ;(*initial temp in the cavity in K*)
icm = m[0] ==0.0045  ;(*initial mass of the vaporized gas. Assuming it just starts at 77k at 1atm and then adding heat*)

sol1 = NDSolve[{TEq,mEq, icT, icm}, {T[t], m[t]},{t, 0, tf} ] ;

P1[t_] = m[t]*R*T[t]/V /.sol1 ; (*Plugging back to get shroud pressure as functon of time*)

(*PLOTS FOR EXPECTED CONDITIONS*)
Plot[P1[t], {t, 0, 1}, PlotRange -> {{0,1}, {103000, 100000}}, GridLines ->Automatic,  ImageSize -> "Large", PlotLabels -> Automatic,  AxesLabel -> {Time (s),Pressure (Pa) }]
Plot[{T[t]/.sol1},{t,0,tf},PlotRange -> {{0,tf}, {200, 0}},GridLines ->Automatic, ImageSize->"Large", PlotLabels->Automatic, AxesLabel -> {Time (s),Temperature (K) }]
Plot[{m[t]/.sol1},{t,0,tf},PlotRange -> {{0,tf}, {0.005, 0}},GridLines ->Automatic, ImageSize->"Large",PlotLabels->Automatic,  AxesLabel -> {Time (s),Mass (Kg) }]

(*SOLUTION FOR DIFFERENT HEAT LOADS, Q: *)
mDotEvap2 = QVary/h ;
TEq = T'[t] == 1/(m[t]*Cp)(mDotEvap2*h-mDotOut*Cp*T[t]-QAmb) ; (*Diff Eq for Temperature in the cavity*)
mEq = m'[t] == mDotEvap2 - mDotOut ; (*Conservation of mass*)
sol2 = Table[NDSolve[{TEq,mEq, icT, icm}, {T[t], m[t]},{t, 0, tf} ],{QVary, {30, 100, 500, 1000}}] ;
P2[t_] = m[t]*R*T[t]/V /.sol2 ; (*Plugging back to get shroud pressure as functon of time*)

Plot[P2[t], {t, 0, 120}, PlotRange -> {{0,120}, {107000, 100000}}, GridLines ->Automatic,  ImageSize -> "Large", PlotLabels ->Automatic,  AxesLabel -> {Time (s),Pressure (Pa) }]
Plot[{T[t]/.sol2},{t,0,tf},PlotRange -> {{0,tf}, {200, 0}},GridLines ->Automatic, ImageSize->"Large", PlotLabels->Automatic, AxesLabel -> {Time (s),Temperature (K) }]
Plot[{m[t]/.sol2},{t,0,tf},PlotRange -> {{0,tf}, {0.005, 0}},GridLines ->Automatic, ImageSize->"Large",PlotLabels->Automatic,  AxesLabel -> {Time (s),Mass (Kg) }]
```

Answer 1

Wrap the first argument in Plot with Evaluate and use the option PlotLegends with your choice of legend labels:

Plot[Evaluate[T[t] /. sol2], {t, 0, tf}, 
 PlotRange -> {{0, tf}, {200, 0}}, GridLines -> Automatic, 
 ImageSize -> "Large",  AxesLabel -> {Time (s), Temperature (K)}, 
 PlotLegends -> LineLegend[Automatic, PromptForm["Q", #] & /@ {30, 100, 500, 1000}, 
   LegendLabel -> T[t]]]

Alternatively, use your choice of colors using the option PlotStyle. For example, add the option PlotStyle -> {Red, Green, Blue, Orange} to get:

Plot[Evaluate[m[t] /. sol2], {t, 0, tf}, 
 PlotRange -> {{0, tf}, {0.005, 0}}, GridLines -> Automatic, 
 ImageSize -> "Large", AxesLabel -> {Time (s), Mass (Kg)}, 
 PlotLegends -> LineLegend[Automatic, PromptForm["Q", #] & /@ {30, 100, 500, 1000}, 
   LegendLabel -> m[t]]]

Add the option PlotStyle -> {Cyan, Magenta, Purple, Gray} to get

Do the same for P2[t]:

Plot[Evaluate@P2[t], {t, 0, 120}, 
 PlotRange -> {{0, 120}, {107000, 100000}}, GridLines -> Automatic, 
 ImageSize -> "Large", AxesLabel -> {Time (s), Pressure (Pa)}, 
 PlotLegends -> LineLegend[Automatic, PromptForm["Q", #] & /@ {30, 100, 500, 1000}, 
   LegendLabel -> "P2(t)"]]

Answer 2

An alternative approch is to use ParametricNDSolveValue once (instead of using NDSolve for each value of QVary):

pndsv = ParametricNDSolveValue[{TEq, m'[t] == QVary/h - mDotOut, icT, 
   icm}, {T[t], m[t]}, {t, 0, tf}, {QVary}]; 

Plot[Evaluate[Table[First@pndsv[q], {q, {30, 100, 500, 1000}}]], {t, 0, 120}, 
 PlotRange -> {{0, 120}, All}, GridLines -> Automatic, 
 ImageSize -> "Large", AxesLabel -> {Time (s), Pressure (Pa)}, 
 PlotLegends -> LineLegend[Automatic, PromptForm["Q", #] & /@ {30, 100, 500, 1000}, 
   LegendLabel -> "T(t)"]]

Plot[Evaluate[Table[Last@pndsv[q], {q, {30, 100, 500, 1000}}]], {t, 0, tf}, 
 PlotRange -> {{0, tf}, {0.005, 0}}, GridLines -> Automatic, 
 ImageSize -> "Large", AxesLabel -> {Time (s), Pressure (Pa)}, 
 PlotLegends -> LineLegend[Automatic, PromptForm["Q", #] & /@ {30, 100, 500, 1000}, 
   LegendLabel -> "m(t)"]]

Plot[Evaluate[Table[Last[pndsv[q]] First[pndsv[q]] R/V, {q, {30, 100, 500, 1000}}]], 
 {t, 0, tf}, 
 PlotRange -> {{0, 120}, {107000, 100000}}, 
 GridLines -> Automatic, ImageSize -> "Large", 
 AxesLabel -> {Time (s), Pressure (Pa)}, 
 PlotLegends -> LineLegend[Automatic, PromptForm["Q", #] & /@ {30, 100, 500, 1000}, 
   LegendLabel -> "P2(t)"]]