def METHOD_NAME(self, variables):
    T_av = variables["X-averaged cell temperature [K]"]
    Q_av = variables["X-averaged total heating [W.m-3]"]
    T_amb = variables["Ambient temperature [K]"]
    y = pybamm.standard_spatial_vars.y
    z = pybamm.standard_spatial_vars.z
    # Account for surface area to volume ratio of pouch cell in surface cooling
    # term
    cell_volume = self.param.L * self.param.L_y * self.param.L_z
    yz_surface_area = self.param.L_y * self.param.L_z
    yz_surface_cooling_coefficient = (
        -(self.param.n.h_cc(y, z) + self.param.p.h_cc(y, z))
        * yz_surface_area
        / cell_volume
    )
    # Edge cooling appears as a boundary term, so no need to account for surface
    # area to volume ratio
    edge_cooling_coefficient = -self.param.h_edge(y, z)
    # Governing equations contain:
    #   - source term for y-z surface cooling
    #   - boundary source term of edge cooling
    # Boundary conditions contain:
    #   - Neumann condition for tab cooling
    # Note: pybamm.source() is used to ensure the source term is multiplied by the
    # correct mass matrix when discretised. The first argument is the source term
    # and the second argument is the variable governed by the equation that the
    # source term appears in.
    self.rhs = {
        T_av: (
            pybamm.laplacian(T_av)
            + pybamm.source(Q_av, T_av)
            + pybamm.source(yz_surface_cooling_coefficient * (T_av - T_amb), T_av)
            + pybamm.source(
                edge_cooling_coefficient * (T_av - T_amb), T_av, boundary=True
            )
        )
        / self.param.rho_c_p_eff(T_av)
    }