def METHOD_NAME(self):
    """
    Returns a function handle, which calculates dy/dt = rhs(y, t).
    Notes
    -----
    ``y`` is a composite array that contains ALL the variables that
    the multiphysics algorithm solves for, e.g., if the constituent
    algorithms are ``TransientFickianDiffusion``, and
    ``TransientFourierConduction``, ``y[0:Np-1]`` refers to the
    concentration, and ``[Np:2*Np-1]`` refers to the temperature
    values.
    """
    def ode_func(t, y):
        # Initialize RHS
        rhs = []
        for i, alg in enumerate(self._algs):
            # Get x from y, assume alg.Np is same for all algs
            x = self._get_x0(y, i)  # again use helper function
            # Store x onto algorithm,
            alg.x = x
            # Build A and b
            alg._update_A_and_b()
            A = alg.A.tocsc()
            b = alg.b
            # Retrieve volume
            V = alg.network[alg.settings["pore_volume"]]
            # Calcualte RHS
            rhs_alg = np.hstack(-A.dot(x) + b)/V
            rhs = np.hstack((rhs, rhs_alg))
        return rhs
    return ode_func