Coffee Cooling ode45 (Exercise 6.1)

function coffee_model()
% Model of a ceramic mug of coffee cooling over time
    radius = 0.04; %m
    height = 0.10; %m
    V = pi * radius^2 * height; % cylindrical cup of coffee
    k = 1.5; % W/m/K
    A = height * 2 * pi * radius + pi * radius ^ 2; % sides and bottom
    d = 0.007; % m
    m = V * 1000; % 1000kg/m^3 density of coffee
    c = 4186; % J/kg/K
    h = 100; % W / m^2 / K
    Tinit = 370; % K
    Troom = 290; % K

    function res = coffee_U_change(t, U)
        % Change in the internal energy of the coffee by conductive and
        % convective heat transfers
        % dependent on the temperature of the coffee
        temp_coffee = U / m / c;
        conduction = k * A / d * (Troom - temp_coffee);
        convection = h * pi * radius ^ 2 * (Troom - temp_coffee);
        res = conduction + convection;
    end

    time_init = 0;
    time_final = 30 * 60; % seconds (30 minutes)

    % Plot the change in temperature over time
    dUdt = @coffee_U_change;
    % Use ode45 to evaluate the change in temperature for each point in
    % time to generate a line. Advanced Euler's method without a constant
    % time step
    [time_2, U_2] = ode45(dUdt, [time_init, time_final], m * c * Tinit);
    plot(time_2, U_2 / m / c, 'b*')

    title('Temperature of Coffee in Ceramic Mug Over Time')
    xlabel('Time (seconds)')
    ylabel('Temperature (Kelvin)')
end