But we still have a considerable difference between a numerical and analytical value, which is not good.
I remove torch.sum from my code and compare the numerical and analytical Jacobian:
# --------------------------------------------------------------------------- #
# without torch.sum()
# --------------------------------------------------------------------------- #
x5 = np.sqrt(2) * s * torch.tensor([0]) + mu
# Weights
omega5 = torch.tensor([1])

# AR(1) shock
aplus = torch.empty(x5.shape)
for epsilon_idx, epsilon_plus in enumerate(x5):
    aplus[epsilon_idx] = a**rho * np.exp(epsilon_plus)


def eval_g_wo_sum(x):
    """ The system of non-linear equilibrium conditions
    without torch.sum()
    x[0]: Capital stock in the next period
    """
    assert len(x) == nvar

    # Consumption today
    con = a * k**alpha * ls**(1-alpha) - x[0]

    # Labor supply tomorrow
    # ls_plus = np.empty(x5.shape)
    ls_plus = torch.empty(x5.shape)
    for aplus_idx, aplus_val in enumerate(aplus):
        ls_plus[aplus_idx] = ls_compute(
            k=x[0], A=aplus_val, alpha=alpha, psi=psi, theta=theta)

    # Capital stock tomorrow
    k_plusplus = torch.empty(aplus.shape)
    for aplus_idx, aplus_val in enumerate(aplus):
        # state_plus = torch.tensor([x[0], aplus_val])[None, :].float()
        k_plusplus[aplus_idx] = torch.tensor(k)  # Needs to be modified

    # Consumption tomorrow
    con_plus = aplus * x[0]**alpha * ls_plus**(1-alpha) - k_plusplus

    # ----------------------------------------------------------------------- #
    # Euler equation
    # ----------------------------------------------------------------------- #
    g0 = 1 / con - beta * alpha * (omega5 * (
        1 / con_plus * aplus * x[0]**(alpha-1) * ls_plus**(1-alpha)))

    print(g0)
    return np.array([g0])


def eval_jac_g_wo_sum(x, flag):
    """ Numerical approximation of the Jacobian of the system of
    non-linear equilibrium conditions
    Jacobian is computed by the automatic gradient of PyTorch """
    assert len(x) == nvar

    row_idx = np.empty(nnzj, dtype=int)  # Row index
    col_idx = np.empty(nnzj, dtype=int)  # Column index

    # Jacobian matrix structure
    if flag:
        for i in range(ncon):
            for j in range(nvar):
                row_idx[j + i * nvar] = i
                col_idx[j + i * nvar] = j

        return (row_idx, col_idx)

    else:
        # Automatic gradient by PyTorch
        assert len(x) == nvar
        _eval_jac_g_wo_sum = np.empty(nnzj)
        for i in range(ncon):
            x_tensor = torch.tensor(x, requires_grad=True)
            eval_g_wo_sum_tensor = eval_g_wo_sum(x_tensor)[i]
            grads = torch.autograd.grad(eval_g_wo_sum_tensor, x_tensor,
                                        create_graph=True)
            for j in range(nvar):
                _eval_jac_g_wo_sum[j + i*nvar] = grads[0][j].detach().numpy()
        return _eval_jac_g_wo_sum


# --------------------------------------------------------------------------- #
# Gradient check: finite difference
# --------------------------------------------------------------------------- #
h = 1e-6
grad = (eval_g_wo_sum(np.array([k + h]))
        - eval_g_wo_sum(np.array([k - h]))) / (2*h)

diff = np.abs(grad - eval_jac_g_wo_sum(np.array([k]), False))
print("Now the difference is {}".format(diff))

but there is still a considerable difference.
I should debug but no clue. Could you show me how should I debug and check how things are going in the automatic gradient in PyTorch?
By the way, I also print g0 and it returns
tensor([-0.7946])
tensor([-0.7947])
tensor([-0.7947], grad_fn=<SubBackward0>)

Is it normal?
Thank for your help.