Big primes enter the picture as follows: Compute the characteristic polynomial $P_n$ of your matrix. Factor $P_n$ over your working prime $p$. This gives you the eigenvalues of the involved Jordan-blocs over $\mathbb F_p$. These eigenvalues are elements of field extensions of degree at most $2n$ of your groundfield $\mathbb F_p$. You can thus get large primes dividing the order of $M$ if they are divisors of $p^k-1$ for $k$ the degree of (at least) one irreducible polynomial (over $\mathbb F_p$) dividing $P(n)$.

A fast way for computing the order of $M(n)$ is thus to compute the characteristic polynomial $P(n)$, factor it over $\mathbb F_p$ and check then if prime-divisors of $p^k-1$ (for $k$ the degree of an involved irreducible polynomial) divide the order.

The order of $M(n)$ over $\mathbb F_p$ is necessarily a divisor of $U=p^{2n-1}\prod_{k=1}^{2n}(p^k-1)$. One can of course remove a factor $p^k-1$ from this product if $P(n)$ mod $p$ has no irreducible divisor modulo $p$ of degree $k$. One can further reduce the size of $U$ by taking greatest common divisors of all involved factors. The exponent $2n-1$ of $p$ in $U$ can of course be replaced by $\mu-1$ where $\mu$ is largest occuring multiplicity among the irreducible divisors of $P(n) \pmod p$. (There seem indeed to be often multiplicities in $P(n)$: This leads to possibly non-trivial Jordan-blocks in $M(n)$ which multiply the order of $M(n)$ by $p^d$ for $d+1$ the dimension of the largest involved Jordan bloc.)

The computational bottleneck of this approach is the factorization of $p^k-1$. The rest is computationally easy using fast exponentiation. Concretetely, You have to check if $M^{U/q}$ is still the identity for $q$ any prime divisor of $U$. If this is the case, replace $U$ by $U/q$ (and recheck if you can remove one more power of $q$ from $U$). Going through all prime divisors of $U$ gives you at the end the exact order of $M(n)$ over $\mathbb F_p$

Big primes enter the picture as follows: Compute the characteristic polynomial $P_L$ of your square matrix $M$ of size $L$. Factor $P_L$ over your working prime $p$. This gives you the eigenvalues of the involved Jordan-blocs over $\mathbb F_p$. These eigenvalues are elements of field extensions of degree at most $L$ of your groundfield $\mathbb F_p$. You can thus get large primes dividing the order of $M$ if they are divisors of $p^k-1$ for $k$ the degree of (at least) one irreducible polynomial (over $\mathbb F_p$) dividing $P_L$.

A fast way for computing the order of $M$ is thus to compute the characteristic polynomial $P_L$ of $M$, factor it over $\mathbb F_p$ and check then if prime-divisors of $p^k-1$ (for $k$ the degree of an involved irreducible polynomial) divide the order.

The order of $M$ over $\mathbb F_p$ is necessarily a divisor of $U=p^{L-1}\prod_{k=1}^{L}(p^k-1)$. One can of course remove a factor $p^k-1$ from this product if $P_L$ mod $p$ has no irreducible divisor modulo $p$ of degree $k$. One can further reduce the size of $U$ by taking greatest common divisors of all involved factors. The exponent $L-1$ of $p$ in $U$ can of course be replaced by $\mu-1$ where $\mu$ is largest occuring multiplicity among the irreducible divisors of $P_L \pmod p$. (There seem indeed to be often multiplicities in $P_L$: This leads to possibly non-trivial Jordan-blocks in $M$ which contribute a factor $p^d$ to the order of $M$ where $d+1$ is the dimension of the largest involved Jordan bloc.)

The computational bottleneck of this approach is the factorization of $p^k-1$. The rest is computationally easy using fast exponentiation. Concretetely, You have to check if $M^{U/q}$ is still the identity for $q$ any prime divisor of $U$. If this is the case, replace $U$ by $U/q$ (and recheck if you can remove one more power of $q$ from $U$). Going through all prime divisors of $U$ gives you at the end the exact order of $M$ over $\mathbb F_p$.

Experimental observation The characteristic polynomial $P_L$ of the square matrix $M=AB$ of even size $L$ seems to be given by $$P_L=(x^2+4)Q^2_{L/2}$$ where $Q_1=1,\quad Q_2=x+2$ and $Q_n=(x+2)Q_{n-1}-xQ_{n-2}$ for $n\geq 3$. (The polynomials $Q_n$ are 'almost' orthogonal polynomials and seem to satisfy the divisibility property $P_a\vert P_b$ if $a\vert b$.) This formula holds for $L\leq 60$ even.