def hessian(f, delta=DELTA):
    """
    Returns numerical hessian function of given input function
    Input: f, scalar function of one or two variables
           delta(optional), finite difference step
    Output: hessian function object
    """
    def hessian_f(*args, **kwargs):
        if len(args) == 1:
            x, = args
            hessianf_x = (
                f(x+delta) + f(x-delta) - 2*f(x)
                )/delta**2
            return hessianf_x
        elif len(args) == 2:
            x, y = args
            if type(x) in [float, int] and type(y) in [float, int]:
                hess_xx = (
                    f(x + delta, y) + f(x - delta, y) - 2*f(x, y)
                    )/delta**2
                hess_yy = (
                    f(x, y + delta) + f(x, y - delta) - 2*f(x, y)
                    )/delta**2
                hess_xy = (
                    + f(x+delta/2, y+delta/2)
                    + f(x-delta/2, y-delta/2)
                    - f(x+delta/2, y-delta/2)
                    - f(x-delta/2, y+delta/2)
                    )/delta**2
                return hess_xx, hess_xy, hess_yy
    return hessian_f