There are infinite elements of S will satisfy the equation.

Suppose X=[a b; c d] satisfies the equation X^2+X+E=0, whre E is the identity matrix.

The eigenvalues λ of X satiesfy λ^2+ λ+1=0, which can be regarded as the characteristic equation of X.

On the other hand, The characteristic equation of X is det(λE-X)=λ^2-(a+d)+ad-bc=0. Hence λ^2+ λ+1=λ^2-(a+d)+ad-bc. Then a+d=-1 and ad-bc=1.

It is obvious a=-d-1, b=1, c=-d^2-d-1 (where d is an arbitrary intege) satisfy the equations a+d=-1 and ad-bc=1.

There are infinite elements of S will satisfy the equation.

Suppose X=[a b; c d] satisfies the equation X^2+X+E=0, whre E is the identity matrix.

The eigenvalues λ of X satiesfy λ^2+ λ+1=0, which can be regarded as the characteristic equation of X.

On the other hand, The characteristic equation of X is det(λE-X)=λ^2-(a+d)+ad-bc=0. Hence λ^2+ λ+1=λ^2-(a+d)+ad-bc. Then a+d=-1 and ad-bc=1.

It is obvious a=-d-1, b=1, c=-d^2-d-1 (where d is an arbitrary intege) satisfy the equations a+d=-1 and ad-bc=1.

Thank you...@Xiaofeng

Yeah Peter, But I was thinking is there anything which is hidden, when I see this question, especially when the field is R and the matrices is made of integers. I just wanted to cross check, everything. that is why

The companion matrix C = [-1  -1;1  0]  and all of its similar matrices  

P^{-1}CP satisfy the equation.

The specified matrix equation leads to a system of linear equations for the elements a, b, c, d. The solution has the following form (for k=1,2,3,...): X1=[-k b; c k-1], X2=[-k -b; -c k-1], X3=[k-1 b; c -k], X4=[k-1 -b; -c -k], X1T, X2T, X3T and X4T where bc=-k2+k-1 (typically, b=k2-k+1 and c=-1 but there are more options for choosing b, c when k2-k+1 is a composite number, e.g. for k=5).

All matrix solutions are similar to the matrix D=[(-1+i\sqrt(3))/2 0; 0 (-1-i\sqrt(3))/2].