@Carl Love 

1.

Yes, it's curious that RowReduce is faster, because for Determinant, the default method is REF.
BTW, replacing in MDS Determinant by RowReduce, in Maple 2019 I get [each after restart]:

CodeTools:-Usage(MDS2(A80,8)):
memory used=1.02GiB, alloc change=32.00MiB, cpu time=15.96s, real time=15.96s, gc time=265.20ms
CodeTools:-Usage(IsMDS(A80,8)):
memory used=2.65GiB, alloc change=60.50MiB, cpu time=15.27s, real time=15.20s, gc time=795.61ms
(In Maple 2018 the results are the same.)

2.

n <= m <= n^2 - n +1

both limits being attained.

@Carl Love 

But of course, if A has integer[], Mod is not necessary, e.g.
A:= LA:-RandomMatrix(64$2, generator= rand(0..1), datatype=integer[kernelopts(wordsize)/8]):

Merry Christmas !

@Magma 

LinearAlgebra:-Modular:-MatrixPower(2, A, k)

@Magma

It also appears in the definition of B[i,j].
I am curious if it takes <1 min. I do not have a 64 MDS matrix to check.

@Magma 

It seems that you have a 32 bit Maple. Please replace integer[8] with integer[4] (twice)

@Magma 

I think you have copied the code while I was editing. Please copy again.

@Carl Love 

Yes. Actually division ring is enough (commutativity is automatic) by Wedderburn's theorem.

@Magma 

Thanks, everything is clear now. So, A is not an arbitrary binary matrix, it is probably derived from a matrix over GF.
I wonder whether this approach is faster than computing the determinants over GF(2^r).
I'll try to find an enhanced version soon.

@Magma

The MDS definition is not in terms of GF. Actually GF is not mentioned at all in the first paper.
That is why I have asked. So, do you know the answers (mainly the first one)?

@Magma 

I did not know about this technique. So, it seems that each r x r block of B corresponds to an element of GF(2^r).

Is it known an embedding of GF(2^r) into the M_r(Z2) ? Or, here it is used another definition for MDS?

@Jjjones98 

I don't see how such a solution could be useful. It is inherently huge (almost each coefficient being a symbolic expression).

Not to mention that for some values of the parameters the system will be inconsistent or undetermined, i.e. the obtained solution is generic only. A complete solution (with all possibilities) would be much much longer.
 

@CyberRob 

f is already in the form you describe. Or, use
collect(f, [k12, k21, Ve], distributed);
It would be better to present mathematically the entire problem (preferably with generic notations: x,y,...  the unknowns,  a,b,... the parameters).

@Carl Love 

OK, but we cannot assume that x-1 and x-3 are the only OP's polynomials.

@Carl Love 

Compare:

(x^2+12)/(x+1) mod 13;

Normal((x^2+12)/(x+1)) mod 13;

@Carl Love 

Probably OP wants

Normal( (x-1)/(x-3) ) mod 13; 
    (x+12)/(x+10)

(the field of fractions).

Or, maybe (in this case)

r := Rem(x-1,x-3,x,'q') mod 13:
q + r/(x-3 mod 13);
    1 + 2/(x+10)