Piecewise function plot with solution from NDSolve

Question

I have these differential equations and I wanted to solve its solutions then plug it into an expression for temperature (which I have listed in the plot below), but the expression for temperature is a piece-wise function one for $t<t_b$ and the other is $t>t_b$.

Mp = 2.4353*10^18;
m = 1.8*10^13;
Rm = (1.8*10^13)/(2.4353*10^18);
v = 13.2 Mp;
Rv = 13.2;
\[Lambda] = 10^-13;
Tb = 10 (m/Mp)^4;
tf = 10^7; 
sol = NDSolve[{dy'[t] + 3 H[t] (1 + 10) dy[t] + \[Lambda] (y[t]^2 - Rv^2) y[t] == 0, H[t] == Sqrt[(0.5 dy[t]^2 + 0.25 \[Lambda] (y[t]^2 - Rv^2)^2 + \[Rho][t])/3], \[Rho]'[t] + 4 H[t] \[Rho][t] == 3 H[t] 10 dy[t]^2, y'[t] == dy[t], y[0] == -30, dy[0] == 9*10^-6, \[Rho][0] == Rm^4}, {y, dy, H, \[Rho]}, {t, 0, tfin},MaxStepFraction -> 0.001, MaxSteps -> 5*10^7]
ysol = y[10^7] /. sol;
ypsol = dy[10^7] /. sol;
Hsol = H[10^7] /. sol
\[Rho][10^7] /. sol;
tb = ((3 Hsol (1 + 10))/(2 \[Lambda] Rv^2)) (Log[(40 Pi^2 Hsol^2 (1 + 10)^2)/(\[Lambda]^2 Rv^6 10)])
Plot[Piecewise[{{Tb Exp[Hsol (x - tb)],x < tb}, {Tb Exp[(\[Lambda] Rv^2 (x-tb))/(6 Hsol (1 + 10))], x > tb}}], {x, 0, 10^7}, PlotRange -> All]

It gave me an error that tb is not a machine sized number, etc. Can anyone help me with this situation? Also, I'm fairly new with mathematica so I've got few more problems, does my code look correct? For example, I substituted "sol" to a function of x, but every quantity needs to be evaluated at some t, so I'm not sure. Just a note, the number 10 you see is something that I can change later so I just left it like that, example (1+10), I did not add it so that it is easier to change later.

Basically, I need to plot $T$ vs. $t$

ADDITIONAL QUESTION: Suppose I have these codes, can anyone help me determine the proper $tf$ and $y[0]$ to get a plot like below, note that in my code $T = \frac{3 \rho[t]}{10 \pi^2}$ where $T$ is the temperature. Also, dy[0]=b shouldn't be too sensitive, $10^{-6}$ would be ok, it is just $y[0]$ and $tf$ that I somehow have some difficulties, although $0 \lt y[0] \lt 50$ (maybe it can go beyond 50 but 500 is too much I think). For some combination of $y[0]$ and $tf$, I'm getting an error saying ""For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions". I'm not sure how to determine the correct combination to get the plot, it's not necessary to get the exact same plot, but $T$ should somehow cut the line $H$ and go below it, like a slanted V.

Mp = 2.4353*10^18;
m = 1.8*10^13;
Rm = (1.8*10^13 )/(2.4353*10^18);
\[Sigma] = 2.24*10^19;
\[Sigma]m = (2.24*10^19)/(2.4353*10^18);
\[Lambda] = ((1.8*10^13)/(2.4353*10^18))^4;
tf = 10^9;
sol[a_, b_] := NDSolve[{dy'[t] + 303 H[t] dy[t] + 4 \[Lambda] (y[t]^2 - \[Sigma]m^2) y[t] == 0, H[t] == Sqrt[(0.5 dy[t]^2 + \[Lambda] (y[t]^2 - \[Sigma]m^2)^2 + \[Rho][t])/3], \[Rho]'[t] + 4 H[t] \[Rho][t] == 300 H[t] dy[t]^2, y'[t] == dy[t], y[0] == -a, dy[0] == b, \[Rho][0] == Rm^4}, {y, dy, H, \[Rho]}, {t, 0, tf}]
sol[1, 1/100000]
Manipulate[Plot[{Evaluate[H[t] /. First@sol[a, b]], Evaluate[(3 \[Rho][t]/10 Pi^2)^(1/4) /. First@sol[a, b]]}, {t, 0, tf}, PlotRange -> Automatic], {{a, 1}, 0, 50, Appearance -> "Labeled"}, {{b, 1/100000}, 0, 0.00009, Appearance -> "Labeled"}]

---

I have resolved the previous problem but suppose I wanted to plot $H$ and $T = (3\rho(t)/10 \pi^2)^{1/4}$ vs. $N$. I have written the code below but when I plot them it just shows the "$N$" axis, no axis for $H$ and $T$. What seems to be wrong in the code? I tried figuring out if I have some syntax error etc. but so far, I don't see anything wrong.

Mp = 2.4353 10^18; 
m = 1.8 10^15; 
Rm = (1.8 10^15 )/(2.4353 10^13); 
\[Sigma] = 2.24 10^19; 
R\[Sigma] = (2.24 10^19)/(2.4353 10^18); 
\[Lambda] = ((1.8 10^15)/(2.4353 10^18))^4; 
tf = 10^8; 
H[t_] = Sqrt[1/3 (0.5 y'[t]^2 + \[Lambda] (y[t]^2 - R\[Sigma]^2)^2 + \[Rho][t])]; 

sol = ParametricNDSolveValue[Rationalize[{y''[t] + 3 (1 + 0.5) H[t] y'[t] + 4 \[Lambda] (y[t]^2 - R\[Sigma]^2) y[t] == 0, \[Rho]'[t] + 4 H[t] \[Rho][t] == 1.5 H[t] y'[t]^2, D[Nx[t], t] == H[t], y[0] == -a, y'[0] == b, \[Rho][0] == Rm^4, Nx[0] == 0}, 0], {y, \[Rho], Nx}, {t, 0, tf}, {a, b}, AccuracyGoal -> 16, MaxSteps -> Infinity] 

Manipulate[ParametricPlot[Evaluate[{{Nx[t], ((3 \[Rho][t])/(10 \[Pi]^2))^(1/4)}, {Nx[t], H[t]}} /. Thread[{y, \[Rho], Nx} -> sol[a, b]]], {t, 0, tf}, PlotRange -> Automatic], {{a, 1}, 0, 2, Appearance -> "Labeled"}, {{b, 1/100000}, 0, 0.00009, Appearance -> "Labeled"}]

Answer 1

Not a answer, but too long for a comment.

Well, as said above, I'm afraid it's hard to determine a and b just from a coding perspective. Put even more bluntly, currently we don't even know whether your modeling is correct or not. Anyway, your coding can be improved, here are some suggestions:

1. NDSolve is able to handle high-order derivative without any difficulty, there's no need to rewrite y''[t] to dy'[t], and as far as I can tell, rewriting y'[t] as dy[t] won't make NDSolve work any better.

2. On the other hand, eliminating H[t] from the equation set will improve the performance of NDSolve, because there's no derivative for H[t] in the equation set, NDSolve will treat the equations as a DAE system. (Indeed, in some cases NDSolve will transform the DAE system to an ODE system, but not always, and the warning about IDA you encountered is exactly an evidence for the triggering of DAE solver, just search IDA in the document. ) The DAE solver of NDSolve is weaker compared to the ODE solver, so it's good to take the expression of H[t] out of NDSolve.

3. You can make use of ParametricNDSolveValue since a and b is under-determined.

4. Since the parameters are extreme, we need to raise the WorkingPrecision of NDSolve. Also, we need to use Rationalize to raise the precision of the parameters because you've used MachinePrecision numbers. (Actually in this case we just need to raise AccuracyGoal to e.g. AccuracyGoal -> 16, see this post for more information. Given AccuracyGoal is a subtle option, I'll stick to WorkingPrecision in the rest part of this post. )

5. The graph you posted is a logarithmic graph, so we can use LogPlot.

6. The position of Evaluate is wrong, it should be outside of the first argument of Plot, or it won't "break through" the HoldAll attribute of Plot.

7. With all due respect, you should put a bit more effort into learning the core language of Mathematica.

The following is an improved version of the code:

Mp = 2.4353 10^18;
m = 1.8 10^13;
Rm = (1.8 10^13)/(2.4353 10^18);
σ = 2.24 10^19;
σm = (2.24 10^19)/(2.4353 10^18);
λ = ((1.8 10^13)/(2.4353 10^18))^4;
tf = 10^9;
H[t_] = Sqrt[1/3 (0.5 y'[t]^2 + λ (y[t]^2 - σm^2)^2 + ρ[t])]
sol = ParametricNDSolveValue[
  Rationalize[{y''[t] + 303 H[t] y'[t] + 4 λ (y[t]^2 - σm^2) y[t] == 0, 
    ρ'[t] + 4 H[t] ρ[t] == 300 H[t] y'[t]^2, y[0] == -a, 
    y'[0] == b, ρ[0] == Rm^4}, 0], {y, ρ}, {t, 0, tf}, {a, b}, 
  WorkingPrecision -> 32]

LogPlot[{H[t], (3/10 ρ[t] π^2)^(1/4)} /. 
   Thread[{y, ρ} -> sol[200, 10^-1]] // Evaluate, {t, 0, tf/100}, PlotRange -> All]

As we can see, under the parameters chosen above, though not crossing curve of $H$, the curve of $T$ bends. Though the goal isn't achieved, it's a small step forward, I think.