Giant Babylonian Matrices
We produce such a matrix of any dimension n. It is simply a matrix of all 1's except 2's on the diagonal. 

r = range(n)
M = [[1 + (i==j) for i in r] for j in r]


Example:

n=5
M=[[2 1 1 1 1]
   [1 2 1 1 1]
   [1 1 2 1 1]
   [1 1 1 2 1]
   [1 1 1 1 2]]

 #Within 10 steps, the procedure converges to the matrix's square root
  [[ 1.28989795  0.28989795  0.28989795  0.28989795  0.28989795]
   [ 0.28989795  1.28989795  0.28989795  0.28989795  0.28989795]
   [ 0.28989795  0.28989795  1.28989795  0.28989795  0.28989795]
   [ 0.28989795  0.28989795  0.28989795  1.28989795  0.28989795]
   [ 0.28989795  0.28989795  0.28989795  0.28989795  1.28989795]]


In general, the procedure converges for matrices with positive real eigenvalues. This is because if we work in the basis of eigenvectors of the original matrix, it and all its iterates are diagonal matrices, and we can view the procedure as applying the Babylonian method to each individual diagonal entry. The starting matrix is the still the identity in this basis.

Since the Babylonian method converges only if the original value was positive, the matrix version converges when every eigenvalue is real and positive. I'm not sure what happens for complex eigenvalues.

I constructed the given matrix to have positive eigenvalues as follows. First, start with the matrix of all ones, which is rank-1, and so has a single eigenvalue (of n) and all others 0. Then, add the identity, which shifts the eigenvalues up by 1 to get (n+1,1,1,...,1). Iterating makes the eigenvalues converge to (sqrt(n+1),1,1,...,1). Note that the 1's are unchanged by each iteration.