Solving multiple differential equations parametrically

Question

I have the following code, where I am trying to solve one set of differential equations parametrically. I then want to use this solution to solve another equation parametrically. However, I get the error: ParametricNDSolveValue::pdvar: Dependent variables {Rast,s} cannot depend on parameters {W,Rast,[Rho],s0,a0}.

The code is as follows:

Clear["Global`*"]

SetOptions[Plot, Frame -> True, Axes -> False, LabelStyle -> {FontFamily -> "Baskerville", FontSize -> 12, Bold}, FrameStyle -> Directive[Thick, Black], 
   ImageSize -> Medium]; 

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]); 

scrit[dens_] := (7.48*Sqrt[dens/2000])/10^4; 

Yarkovsky

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]*Qpr)/(c^2*(32*Pi*Rast*\[Rho]*a[t]^2*Sqrt[1 - e[t]^2]))); 

YdadtR\[Rho]a = -((1/c)*((3*L[t])/(16*Pi*Rast*\[Rho]*(1 - e[t]^2)*Sqrt[G*Mwd*a[t]]))); 

YdedtR\[Rho]a = -((1/c)*((3*L[t]*(1 - Sqrt[1 - e[t]^2]))/(32*Pi*Rast*\[Rho]*e[t]*Sqrt[G*Mwd*a[t]^3]))); 

dadtR\[Rho]a = Evaluate[YdadtR\[Rho]a*UnitStep[Rast - 1] + RDdadtR\[Rho]a]; 

dedtR\[Rho]a = Evaluate[YdedtR\[Rho]a*UnitStep[Rast - 1] + RDdedtR\[Rho]a]; 

YORPdsdtWR\[Rho]s = (W/(2*Pi*\[Rho]*Rast^2))*(1/(a[t]^2*Sqrt[1 - e[t]^2]))*(U*(L[t]/Lsun)); 

solR\[Rho]a = ParametricNDSolveValue[{Derivative[1][s][t] == YORPdsdtWR\[Rho]s, Derivative[1][a][t] == dadtR\[Rho]a, Derivative[1][e][t] == dedtR\[Rho]a, s[0] == s0, a[0] == a0, 
     e[0] == 0.3}, {s, a, e}, {t, 0, 9*Gyr}, {W, Rast, \[Rho], s0, a0}]; 

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

Manipulate[Column[{Style["YORP working plot", Bold], Plot[fun[func, t], {t, 0, 9*Gyr}, FrameTicks -> fticks, 
     Epilog -> {Red, Dashed, InfiniteLine[{{0, scrit[\[Rho]]}, {10, scrit[\[Rho]]}}]}, PlotStyle -> {Directive[Blue, Thickness[0.01]]}], Style["Compiled Plot", Bold], 
    If[comp === {}, Plot[fun[func, t], {t, 0, 9*Gyr}, FrameTicks -> fticks, Epilog -> {Red, Dashed, InfiniteLine[{{0, scrit[\[Rho]]}, {10, scrit[\[Rho]]}}]}, 
      PlotStyle -> {Directive[Blue, Thickness[0.01]]}], Plot[comp, {t, 0, 9*Gyr}, FrameTicks -> fticks, 
      Epilog -> {Red, Dashed, InfiniteLine[{{0, scrit[\[Rho]]}, {10, scrit[\[Rho]]}}]}, PlotStyle -> {Directive[Blue, Thickness[0.01]]}]]}], {{func, 1}, {1 -> "s"}}, 
  {{W, 0}, 0, 1, 0.1, Appearance -> "Labeled"}, {{Rast, 50}, 1, 4000, 10, Appearance -> "Labeled"}, {{\[Rho], 2000}, 1000, 7000, 50, Appearance -> "Labeled"}, 
  {{s0, 0, "s0 (\!\(\*SuperscriptBox[\(rads\), \(-1\)]\))"}, 0, 100, 1, Appearance -> "Labeled"}, {{a0, 1, "a0 (au)"}, 1, 20, 1, Appearance -> "Labeled"}, 
  Button["Append", AppendTo[comp, fun[func, t]]], Button["Reset", comp = {}], TrackedSymbols :> {func, W, Rast, \[Rho], s0, a0}, 
  Initialization :> {comp = {}, fun[func_, t_] := solR\[Rho]a[W, Rast, \[Rho], s0, a0][[1]][t]}]

Any help would be appreciated.

Answer 1

The code has three difficulties. First, asol appears in YORPdsdtWRρs as asol[Rast, ρ, a0, t], which evaluates to

Clear[asol]
asol[Rast_, ρ_, a0_, t_] := solRρa[Rast, ρ, a0][[1]][t]
asol[Rast, ρ, a0, t]
(* Rast[t] *)

To prevent this, instead use

Clear[asol]
asol[Rast_?NumericQ, ρ_?NumericQ, a0_?NumericQ, t_] := solRρa[Rast, ρ, a0][[1]][t]
asol[Rast, ρ, a0, t]
(* asol[Rast, ρ, a0, t] *)

which does not evaluate prematurely. Next, plot asol with the default parameter values given in Manipulate in the queston.

Plot[Evaluate[asol[50, 2000, 1, t]], {t, 0, 1.726 10^-6}, 
    ImageSize -> Large, AxesLabel -> {t, a}, LabelStyle -> {15, Bold, Black}]

Visibly, the stellar radius goes to zero at about t = 1.7260 10^-6, which is the source of the error message from Evaluate[asol[50, 2000, 1, t]].

Finally, ParametricNDSolve objects with attempting to evaluate solWRρs, because "Dependent variables {s,asol[Rast,ρ,a0,t]} cannot depend on parameters {W,Rast,ρ,s0,a0}". To solve this final problem, merge the computations of solRρa and solWRρs.

YORPdsdtWRρs = (W/(2*Pi*ρ*Rast^2))*(1/(a[t]^2*Sqrt[1 - ecc^2]))*(U*(L[t]/Lsun));
solWRρs = ParametricNDSolveValue[{Derivative[1][s][t] == YORPdsdtWRρs, s[0] == s0, 
    Derivative[1][a][t] == dadtRρa, Derivative[1][e][t] == dedtRρa, a[0] == a0, 
    e[0] == 0.3}, {s[t], a[t], e[t]}, {t, 0, 9*Gyr}, {W, Rast, ρ, s0, a0}];

(Note that I have not corrected the second problem by changing the upper bound on t immediately above, because it causes no substantial harm.)

Plot[Evaluate[First@solWRρs[1, 50, 2000, 0, 1]], {t, 0, 1.726 10^-6}, 
    ImageSize -> Large, AxesLabel -> {t, s}, LabelStyle -> {15, Bold, Black}]

The plot of a, and of e as well, can be recovered by

Plot[Evaluate[Rest@solWRρs[1, 50, 2000, 0, 1]], {t, 0, 1.726 10^-6}, 
    ImageSize -> Large, AxesLabel -> {t, a}, LabelStyle -> {15, Bold, Black}]