# equation solving – Using Solution to Differential Eqution 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 just got get any output on the pressure plot (the last one).

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