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

I am using [`NDSolve`](https://reference.wolfram.com/language/ref/NDSolve.html) 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`](https://reference.wolfram.com/language/ref/NDSolve.html) 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]
```