Colour each graph differently

Question

I have a Manipulate function, which allows me to compile plots on a single plot (see below). I am trying to make each line on the compiled plot to be a different colour, however, Mathematica does not allow me to do this. What am I doing wrong?

In the code below, I have tried to ensure the first of the compiled plots is blue and the second is red, however, they both turn out red.

(* Constants *)

au = QuantityMagnitude[UnitConvert[Quantity[1, "AstronomicalUnit"], "Meters"]]; 
c = QuantityMagnitude[UnitConvert[Quantity[1, "SpeedOfLight"], "MetersPerSecond"]]; 
Qpr = 1; 
Lsun = QuantityMagnitude[UnitConvert[Quantity[1, "SolarLuminosity"], "Watts"]]; 
Rsun = QuantityMagnitude[UnitConvert[Quantity[1, "SolarRadius"], "Meters"]]; 
Msun = QuantityMagnitude[UnitConvert[Quantity[1, "SolarMass"], "Kilograms"]]; 
G = QuantityMagnitude[UnitConvert[Quantity[1, "GravitationalConstant"], ("Meters"^2*"Newtons")/"Kilograms"^2]]; 
year = QuantityMagnitude[UnitConvert[Quantity[1, "Years"], "Seconds"]]; 
Myr = year*10^6; 
Gyr = year*10^9; 
Mwd = 0.6*Msun; 
Cst = 1.27; 
U = 1*10^17; 

(* Functions *)

L[t_] := (3.26*Lsun*(Mwd/(0.6*Msun)))/(0.1 + t/Myr)^1.18; 
Roche[dens_] := (0.65*Cst*Rsun*(Mwd/(0.6*Msun))^(1/3))/(dens/3000)^3^(-1); 
Papsis[t_] := a[t]*(1 - e[t]); 

(* Radiative Drag *)

RDdadtRρa = -((3*L[t]*Qpr*(2 + 3*e[t]^2))/(c^2*(16*Pi*ρ*Rast*a[t]*(1 - e[t]^2)^(3/2)))); 
RDdedtRρa = -((15*L[t]*e[t])/(c^2*(32*Pi*Rast*ρ*a[t]^2*Sqrt[1 - e[t]^2]))); 

RDsolRρa = 
  ParametricNDSolveValue[
    {Derivative[1][a][t] == RDdadtRρa, 
     Derivative[1][e][t] == RDdedtRρa, 
     a[0] == a0, e[0] == 0.3
    }, {a, e}, {t, 0, 9*Gyr}, {Rast, ρ, a0}]; 

fRDticks = {{Automatic, Automatic}, {Charting`FindTicks[{0, 1}, {0, 1/Myr}], Automatic}}; 

Manipulate[
  Column[{
    Style["Radiative Drag Working Plot", Bold], 
    Plot[
      fun[func, t]/scale[func], {t, 0, 9*Gyr}, 
      FrameTicks -> fRDticks, 
      Epilog -> {Red, Dashed, InfiniteLine[{{0, Roche[ρ]}, {10, Roche[ρ]}}]}, 
      PlotStyle -> {Directive[Blue, Thickness[0.01]]}
    ],
    Style["Compiled Plot", Bold], 
    If[comp === {}, 
      Plot[fun[func, t]/scale[func], {t, 0, 9*Gyr}, FrameTicks -> fRDticks, Epilog -> {Red, Dashed, InfiniteLine[{{0, Roche[ρ]}, {10, Roche[ρ]}}]}, PlotStyle -> {Directive[Blue, Thickness[0.01]]}], 
      Plot[comp/scale[func], {t, 0, 9*Gyr}, FrameTicks -> fRDticks, Epilog -> {Red, Dashed, InfiniteLine[{{0, Roche[ρ]}, {10, Roche[ρ]}}]}, PlotStyle -> {Directive[Blue, Thickness[0.01]], Directive[Red, Thickness[0.01]]}]
    ]
  }], 
  {{func, 1}, {1 -> "a", 2 -> "e", 3 -> "q"}}, 
  {{Rast, 0.005}, 0.0001, 0.1, 0.001, Appearance -> "Labeled"}, 
  {{ρ, 3000}, 1000, 7000, 50, Appearance -> "Labeled"}, 
  {{a0, 10, "a0 (au)"}, 1, 20, 0.2, Appearance -> "Labeled"}, 
  Button["Append", AppendTo[comp, fun[func, t]]], 
  Button["Reset", comp = {}], 
  TrackedSymbols -> {func, Rast, ρ, a0}, 
  Initialization :> {
     comp = {}, 
     fun[sel_, t_] := Switch[sel, 1, RDsolRρa[Rast, ρ, a0*au][[1]][t], 2, RDsolRρa[Rast, ρ, a0*au][[2]][t], 3, RDsolRρa[Rast, ρ, a0*au][[1]][t]*(1 - RDsolRρa[Rast, ρ, a0*au][[2]][t])], 
     scale[sel_] := Switch[sel, 1 | 3, au, 2, 1]
  }
]

Answer 1

You need to set PlotStyle -> {Directive[Red, Thickness[0.01]]} for the second plot. Here is the changes code:

fRDticks = {{Automatic, 
    Automatic}, {Charting`FindTicks[{0, 1}, {0, 1/Myr}], Automatic}};

Manipulate[
 Column[{Style["Radiative Drag Working Plot", Bold], 
   Plot[fun[func, t]/scale[func], {t, 0, 9*Gyr}, 
    FrameTicks -> fRDticks, 
    Epilog -> {Red, Dashed, 
      InfiniteLine[{{0, Roche[\[Rho]]}, {10, Roche[\[Rho]]}}]}, 
    PlotStyle -> {Directive[Blue, Thickness[0.01]]}],
   Style["Compiled Plot", Bold],
   If[comp === {},
    
    Plot[fun[func, t]/scale[func], {t, 0, 9*Gyr}, 
     FrameTicks -> fRDticks, 
     Epilog -> {Red, Dashed, 
       InfiniteLine[{{0, Roche[\[Rho]]}, {10, Roche[\[Rho]]}}]}, 
     PlotStyle -> {Directive[Red, Thickness[0.01]]}],
    
    Plot[comp/scale[func], {t, 0, 9*Gyr}, FrameTicks -> fRDticks, 
     Epilog -> {Red, Dashed, 
       InfiniteLine[{{0, Roche[\[Rho]]}, {10, Roche[\[Rho]]}}]}, 
     PlotStyle -> {Directive[Red, Thickness[0.01]], 
       Directive[Red, Thickness[0.01]]}]
    ]}], {{func, 1}, {1 -> "a", 2 -> "e", 3 -> "q"}}, {{Rast, 0.005}, 
  0.0001, 0.1, 0.001, Appearance -> "Labeled"}, {{\[Rho], 3000}, 1000,
   7000, 50, Appearance -> "Labeled"}, {{a0, 10, "a0 (au)"}, 1, 20, 
  0.2, Appearance -> "Labeled"}, 
 Button["Append", AppendTo[comp, fun[func, t]]], 
 Button["Reset", comp = {}], 
 TrackedSymbols -> {func, Rast, \[Rho], a0}, 
 Initialization :> {comp = {}, 
   fun[sel_, t_] := 
    Switch[sel, 1, RDsolR\[Rho]a[Rast, \[Rho], a0*au][[1]][t], 2, 
     RDsolR\[Rho]a[Rast, \[Rho], a0*au][[2]][t], 3, 
     RDsolR\[Rho]a[Rast, \[Rho], a0*au][[1]][
       t]*(1 - RDsolR\[Rho]a[Rast, \[Rho], a0*au][[2]][t])], 
   scale[sel_] := Switch[sel, 1 | 3, au, 2, 1]}]

Answer 2

I'd recommend avoiding boilerplate repetitions. Since your plotting options for the two branches of your If are the same, and only the plotting functions change, then instead of If[condition, Plot[], Plot[]], I would recommend Plot[If[condition, expression1, expression2], repeataedIdenticalOptions].

In other words, change your If[....] to:

Plot[
  If[comp === {}, fun[func, t]/scale[func], comp/scale[func]], 
  {t, 0, 9*Gyr}, 
  FrameTicks -> fRDticks, 
  Epilog -> {Red, Dashed, InfiniteLine[{{0, Roche[ρ]}, {10, Roche[ρ]}}]}, 
  PlotStyle -> Directive[Red, Thickness[0.01]]
]