def clmixhess(obj, exe, arg1, arg2, delta=DELTA):
    """
    Returns numerical mixed Hessian function of given class method
    with respect to two class attributes
    Input: obj, general object
           exe (str), name of object method
           arg1(str), name of object attribute
           arg2(str), name of object attribute
           delta(float, optional), finite difference step
    Output: Hessian function object
    """
    f, x = get_method_and_copy_of_attribute(obj, exe, arg1)
    _, y = get_method_and_copy_of_attribute(obj, exe, arg2)
    def hess_f(*args, **kwargs):
        hess_val = numpy.zeros(x.shape + y.shape)
        it = numpy.nditer(x, op_flags=['readwrite'], flags=['multi_index'])
        for xi in it:
            i = it.multi_index
            jt = numpy.nditer(y, op_flags=['readwrite'], flags=['multi_index'])
            for yj in jt:
                j = jt.multi_index
                xi += delta/2
                yj += delta/2
                fpp = f(*args, **kwargs)
                yj -= delta
                fpm = f(*args, **kwargs)
                xi -= delta
                fmm = f(*args, **kwargs)
                yj += delta
                fmp = f(*args, **kwargs)
                xi += delta/2
                yj -= delta/2
                hess_val[i + j] = (fpp + fmm - fpm - fmp)/delta**2
        return hess_val
    return hess_f