def jacobi(A, b, x0, maxiter=10**3, tol=1e-6):
    x = x0.copy()
    diagA = np.diag(A)
    d = b - A @ x
    k = 0

    while True:
        x = x + np.divide(d, diagA)
        d = b - A @ x
        k += 1

        if np.linalg.norm(d) < tol or k > maxiter:
            break

    return x, k