Using Solution to Differential Equation in New Function

I am trying to feed the solutions from my differential equations into another function.

I am using NDSolve to solve two differential equations simultaneously. It seems to be working ok based on the plots it is outputting (I need to work on the physical interpretation, but at this point, I am just trying to get this to output pressure as a function of time on the plot).

Now I want to take the result of the NDSolve and plug it back into the pressure equation to see how pressure changes over time. I assumed I could just put in the interpolated numerical solution directly into the equation and then plot it. I have tried a lot of different stuff but I am new to Mathematica and am probably doing a bunch of stuff wrong.

I am not getting any errors but I am just not getting any output on the pressure plot (the last plot in the code).

Note that the pressure equation P2 (the one I am trying to evaluate using diff eq sol’s) is actually used as P to solve the differential equations in the first place (i.e. P2 = P).

Here is where I am currently at with the code:

(* This calculation is intended to check the pressure build up inside the LN2 cavity*)

ClearAll["Global`*"]

Q =  100 ;(*Heat into the shroud W/m^2*)
QAmb = 0  ;(*Heat loss to ambinet*)
A =  1.935*10^-5 ;(*Area of orifice m^2  Based on 1/4 inch pipe with 0.028 inch wall thickness*)
h =  199 ;(*Heat of vaporization of LN2 kJ/kg*)
Cp =  1.039 ;(*Specific heat of LN2 kJ/(kg K) *)
R =  0.2968 ;(*Gas constant for nitrogen kJ/(kg K)*)
γ = 1.40  ;(*Specific heat ratio*)
V = 0.001 ; (*Enclosed volume m^3*)
icT = T[0] == 300 ;(*initial temp in the cavity*)
icT2 = T'[3600] == 0 ;(*After a long time we will reach steady state and T will stop changing*)
icm = m[0] == 0.001 ;(*initial mass of the vaporized gas*)
tf = 10 ;(*Final time in seconds*)

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

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*)

sol1 = NDSolve[{TEq,mEq, icT, icm}, {T[t], m[t]},{t, 0, tf} ] ;
 P2[t] = ( m[t]/.sol1)*R*(T[t]/.sol1)/V  ; (*Plugging back to get shroud pressure as functon of time*)


Plot[{T[t]/.sol1},{t,0,tf}, 
PlotRange -> {{0,tf}, {800, 0}}, 
ImageSize->"Large",
PlotLabels->Automatic]

Plot[{m[t]/.sol1},{t,0,tf}, 
PlotRange -> {{0,tf}, {1, 0}}, 
ImageSize->"Large",
PlotLabels->Automatic]

Plot[P2[t],{t,0,tf}, 
PlotRange -> Automatic, 
ImageSize->"Large",
PlotLabels->Automatic]

P2[t_] = m[t]*R*T[t]/V /. sol1;
Plot[P2[t], {t, 0, tf}, PlotRange -> Automatic, ImageSize -> "Large", 
 PlotLabels -> Automatic]