Newton's law of cooling

You are supposed to serve a cup of hot coffee and waiting eagerly for getting it little cold so that the most desired sip in the coffee cup can be made. You must be knowing that the equation that governs the cooling down of hot coffee is the Newton's law of cooling, given by $$\frac{dT}{dt}= -k(T-T_s)$$ Here, $T$ is the temperature of the body and $T_s$ is the temperature of the surrounding. $k$ is a constant of proportionality, called heat transfer coefficient.

This is a first order ODE, which can be solved by any standard procedure. Let us try to solve the problem by odieint function from SciPy package, and also by RK4 method.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# basic parameters
T0 = 70
Ts = 30
k_list  = np.array([0.005, 0.05, 0.01])
color   = ['red', 'blue', 'green']

# Discretize the independent variable 't'
t_init, t_final, step_size = 0, 300, 0.1
t = np.arange(t_init, t_final, step_size)

@np.vectorize
def model(T, t, k):
    '''
    dT/dt = - k(T-Ts)
    '''
    slope = - k*(T-Ts)
    return slope

# solving by 'odeint' function
for i in range(len(k_list)):
    sol = odeint(model, T0, t, args=(k_list[i], ))
    T = sol[:, 0]
    plt.plot(t, T, color=color[i], label='k = %0.3f'%k_list[i])
plt.xlabel('$t$', fontsize=14)
plt.ylabel('$\\frac{dT}{dt}$', fontsize=14)
plt.legend()
plt.show()

# Create array for dependent variable "T"
T = np.zeros(len(t))

T[0] = T0

# Defining the helper functions
    
def rk4(T, t, h, k):
    k1 = model(T, t, k)
    k2 = model(T+k1*h/2, t+h/2, k)
    k3 = model(T+k2*h/2, t+h/2, k)
    k4 = model(T+k3*h, t+h, k)
    k_avg = (k1+2*(k2+k3)+k4)/6
    T_new = T + k_avg*h
    return T_new

h = step_size
for j in range(len(k_list)):
    for i in range(len(t)-1):
        T[i+1]=rk4(T[i], t[i], h, k_list[j])
    plt.plot(t, T, color=color[j], label='k = %0.3f'%k_list[j])
    plt.xlabel('$t$', fontsize=14)
    plt.ylabel('$\\frac{dT}{dt}$', fontsize=14)

for i in range(len(k_list)):
    plt.text(150, 30+i*30/3, '$k_{%d} = $'%(i+1), fontsize=13)
    
for i in range(len(k_list)):
    plt.text(180, 30+i*30/3, ' %0.3f '%k_list[i], fontsize=13)

plt.show()

So, how quickly a coffee cup will cool down depends on the value of $k$.