plotting – Colour each graph differently

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(Rho)a = -((3*L(t)*Qpr*(2 + 3*e(t)^2))/(c^2*(16*Pi*(Rho)*Rast*a(t)*(1 - e(t)^2)^(3/2)))); 
RDdedtR(Rho)a = -((15*L(t)*e(t))/(c^2*(32*Pi*Rast*(Rho)*a(t)^2*Sqrt(1 - e(t)^2)))); 

RDsolR(Rho)a = ParametricNDSolveValue({Derivative(1)(a)(t) == RDdadtR(Rho)a, Derivative(1)(e)(t) == RDdedtR(Rho)a, a(0) == a0, e(0) == 0.3}, {a, e}, {t, 0, 9*Gyr}, 
    {Rast, (Rho), 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((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(Blue, 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(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"}, {{(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)})

Any help would be appreciated.