@adel-00 The line

DeltaE:=2*(j*(spin*neighbours)+B*spin:

still has unbalanced parentheses. There are several places the closing parenthesis can go. I can't tell where you want it. Also, if you mean standard matrix multiplication, then spin*neighbours should be spin.neighbours.

The line 

p=exp(-DeltaE):

should be

p:= LinearAlgebra:-MatrixExponential(-DeltaE)

(unless you mean for the exp to be applied elementwise).

There are other errors, but I need to go now. Please don't waste my time: I'm sure that you can find unbalanced parentheses on your own.

@adel-00 Now you have unbalanced parentheses in the line that begins DeltaE:=. There are numerous other errors after that. Ordinary Matrix-Matrix multiplication is done with . instead of *. To add all the elements of a Matrix M use (ArrayTools:-AddAlongDimension@@2)(M).

@adel-00 The above code is Doug's, not mine. Like I said before you skipped a line of his code. You need to insert:

spin:= Spin(N);

Or you can use my code.

Regarding seeing an image of the Matrix: Right-click on the Matrix in the worksheet output to bring up the context menu. Then select Browse. Then click on the Image tab.

@testht06 I'm just curious: What application do you have for finding fractional powers of matrices over prime-power fields? Is it coding theory? What will be the size of the matrices and the degree of the field when you're doing practical problems?

@testht06 I don't anticipate the degree of the characteristic polynomial or its splitting field to be any limitation. However, the program is not set up to deal with the following situations:

1. The characteristic polynomial has multiple nonlinear irreducible factors.
2. There are repeated eigenvalues.
3. The nth root of an eigenvalue doesn't exist in the splitting field.

@Preben Alsholm I mean without the extra parentheses around CCnew, as in `<|>`(-CCnew||(0..dmax-1)). It generates a Matrix of Matrices rather than a "flat" Matrix.

I found that by declaring a parameter in the initial isolve, I could avoid the need for subs in the second isolve. I edited the Answer to refect this simplified technique.

@adel-00 You skipped one line of Doug's code:

spin:= Spin(3);

which you'll want to change to

spin:= Spin(N);

@spalinowy To do multiple equations, simply include all the equations and all the initial conditions in the curly braces that delimit the first argument to dsolve. No other changes are needed. The order doesn't matter.

I just computed the fourth root of your second matrix (from your first Reply above):

M:= <x^3+x^4+x^6+x^7+x, 1+x^2+x^4+x^5+x^6, 1+x+x^2+x^3, x^7+x^6+x^5+x^4;
x^7+x^5+x^4+x^3, x^6+x^4+x^2+1, x^3+x^4+x^6+x^7+x^2+x+1, 1+x^2+x^3+x^4+x^5;
x^7+x^5+x^2, x^7+x^5+x^3+x^2+1, x+x^2+x^6, x^2+x^3+x^5;
x^3+x^4+x^6+1, 1+x^2+x^3+x^4, x^6+x^5+x^4+x^3, x^7+x^3 >;

over GF(2) / (x^8+x^5+x^3+x+1). In this case, the characteristic polynomial is irreducible, so the splitting field has degree 32. The computations are nearly instantaneous, so the Groebner:-Basis technique is much, much better than the Primfield technique.

@Doug Meade A Matrix expression such as

< spin[.., -1] | spin[.., 1..-2] >

can be shortened to

spin[.., [-1, 1..-2]]

I don't know about the relative efficiency or about how ArrayTools:-CircularShift compares.

@Zihan You need to freeze the diff first, or else it will turn to 0 when you freeze the x(t).

@Doug Meade Your technique assumes that once fsolve fails for a particular value of T that it'll fail for any smaller T. Although that seems likely to be true, I don't think that it's a wise policy to assume it.

@Christian Wolinski 

If I make the call

Groebner[Basis]([L], plex(X,Q), characteristic= 2);

then I get results identical to yours:

[Q^24+Q^23+Q^22+Q^21+Q^19+Q^16+Q^14+Q^13+Q^12+Q^10+Q^9+Q^8+Q^6+Q^5+Q^3+Q+1, Q^22+Q^21+Q^19+Q^17+Q^16+Q^15+Q^14+Q^9+Q^7+Q^6+Q^5+Q^2+Q+X]

I wasn't aware that the term order could change the final result. I thought that it only affected the speed with which the result was obtained.

@Markiyan Hirnyk 

Your methods work on a polynomial, but not on more complicated expressions such as

sin(1+eps+eps^2+eps^3)