Fråga om differentialekvation - avsvalning av kaffe

Hej,

Jag har tagit fram en differentialekvation för avslaningen av kaffe med och utan lock. Detta med hjälp av ekvationen 

d(mCpT)/dt = -Qvägg - Qstrålning-Qvätska-Qvap.

På python ser det ut såhär med mina värden:

def avsvalningTu(t,T):
    r2 = 0.0375 #m
    h = 0.08 #m
    r1 = 0.0275 #m
    T2 = 25 #°C
    k = 6.5*10**-3 #m/s (massöverföringstal)
    b = 0.001
    alfa1 = 3 #[W/m**2 °C]
    alfa2 = 200 #[W/m**2 °C]¨
    boltzk = 5.65 * 10**-8 #W/m**2k**4
    emissP = 0.93
    lambd = 0.18
    A = 10.19625
    B = 1730.630
    C = 233.426
    deltaHvap = 22060 #J/mol
    R = 8.3145 #J/mol*k
    Mv = 18.01528*10**-3
    m = 0.18 #kg
    Cp = 4.18*1e3 #J/kg*C
   
    dT1dt = -((emissP*boltzk)*2*np.pi*r1*h*((T+273.15)**4-(T2+273.15)**4))/(m*Cp)-(((1/alfa1)+ (1/alfa2))**(-1)*np.pi*r2*(T-T2))/(m*Cp)-(1/alfa1 + b/lambd + 1/alfa2)**(-1)*2*np.pi*r1*h*(T-T2)/(m*Cp) -(deltaHvap*np.pi*(r2**2)*k/R*(10**(A-(B/(T+C)))/(T+273.15) - 0.5*(10**(A-(B/(25+C)))/(25+273.5))))/(m*Cp)

    return dT1dt
 #~~~ Lös den dynamsika modellen ~~~
T0 = [91] # startemperaturen
t_span = [0,90*60]
sol_dyn = solve_ivp(lambda t, T: avsvalningTu(t,T), t_span, T0)
t = sol_dyn.t
T = sol_dyn.y.T
# ~~~ Hitta den stationära punkten ~~~
T_g = 25 # gissninsvärde
sol_stat = root(lambda T: avsvalningTu(0, T), T_g)
T_stat = sol_stat.x
print(T_stat)
plt.plot(t/60, T, '--', label='Simulering - avsvalning kaffe utan lock')
plt.xlabel('Tid [min]')
plt.ylabel(r'Temperatur $[°C]$')

def avsvalningTm(t,T):
    r2 = 0.0375 #m
    h = 0.08 #m
    r1 = 0.0275 #m
    T2 = 25 #°C
    k = 6.5*10**-3 #m/s (massöverföringstal)
    b = 0.001
    alfa1 = 3 #[W/m**2 °C]
    alfa2 = 200 #[W/m**2 °C]¨
    boltzk = 5.65 * 10**-8 #W/m**2k**4
    emissP = 0.93
    lambd = 0.18
    A = 10.19625
    B = 1730.630
    C = 233.426
    R = 8.3145 #J/mol*k
    Mv = 18.01528*10**-3
    m = 0.18 #kg
    Cp = 4.18*1e3 #J/kg*C
   
    dT1dt = -((emissP*boltzk)*2*np.pi*r1*h*((T+273.15)**4-(T2+273.15)**4))/(m*Cp)-(((1/alfa1)+ (1/alfa2))**(-1)*np.pi*r2*(T-T2))/(m*Cp)-(1/alfa1 + b/lambd + 1/alfa2)**(-1)*2*np.pi*r1*h*(T-T2)/(m*Cp)

    return dT1dt

# ~~~ Lös den dynamsika modellen ~~~
T0 = [92] # startemperaturen
t_span = [0,90*60]

dTdt_Tu=avsvalningTu(0,T0[0])
print(dTdt_Tu)

dTdt_Tm=avsvalningTm(0,T0[0])
print(dTdt_Tm)

sol_dyn = solve_ivp(lambda t, T: avsvalningTm(t,T), t_span, T0)
t = sol_dyn.t
T = sol_dyn.y.T
# ~~~ Hitta den stationära punkten ~~~
T_g = 25 # gissninsvärde
fun_stat = lambda T: avsvalningTm(0, T)
sol_stat = root(fun_stat, T_g)
T_stat = sol_stat.x
print(T_stat)
plt.plot(t/60, T, '^:g', label=r'Simulering - avsvalning kaffe med lock')
plt.xlabel('Tid [min]')
plt.ylabel(r'Temperatur $[°C]$')
plt.legend(loc='best')
plt.grid()
plt.show()

Kravet på kaffekoppen är att den ska hålla 50 grader i minst 20 minuter. Jag undrar nu hur man kan göra för att ta reda på hur mycket serveringstemperaturn på kaffet kan sänkas och och fortfarande klara temperaurkravet. Kan det vara så enkelt som att man ska sätta in värden och studera kurvan? Jag tänker att man bör leta efter det T som gör att arean under kurvan (primTstart-prim50) är 20.  Dock känns det väldigt svårt att hitta den primitva funktionen till en funktion som har så många amdra variabler.