Let $M \in F_p^{n\times n}$ be a square matrix of dimension $n\times n$ with entries over $F_p$.
Then we want to obtain $M$ from $M^2=M\cdot M$. Initially I think of the old diagonalisation procedure of solving the square root of a matrix by:
$$M^2 = P\cdot D \cdot P^{-1}$$
$$M = P\cdot D^{\frac{1}{2}} \cdot P^{-1}$$
Where $D$ is a diagonal matrix with the eigenvalues as entries defined over $F_p$. So $2^{-1} \pmod p$ must be a unit on $F_p^{*}$, this is $2$ must be an unit on $F_p^{*}$ so it has an inverse.
Moreover, when taking $F_2$ as the field, the diagonalisation method doesn't yield the square root matrix as $2$ is not an unit on $F_2^{*}$.
I know that if $M^k \in GL(n,F_2)$ then $\gcd(k,\vert GL(n,F_2) \vert)$ must be equal to $1$ so $M^{k\cdot d}=M$ where $d$ is the multiplicative inverse of $k$ modulo the order of the given general linear group. It works also with the multiplicative order of $M$ in $GL$.
But what if $k=2$ and the field is $F_2$? Is there any method for computing the square root matrix of $M^2$ over $F_2$?
Depends.
You already start with a rather strong assumption, that your matrix $M^2$ can be diagonalized. That is wrong in general and only true for some matrices.
But let's assume that you can write $$M^2 = PDP^{-1}.$$
Then every diagonal entry of $D$ is either $0$ or $1$, and both of these have unique square roots. In fact, you get $D^{1/2} = D$ and thus $M = M^2$. That is, of course, a boring case, but once again shows that your assumption that $M$ is diagonizable is already very restrictive; in the case of the binary field, it means that $M^k = M$ for all powers $k$.
In general it seems you are slightly confused with some things. To be able to get a unique $k$-th root, you need that the map $$f : \mathbb{F}_p \to \mathbb{F}_p, x \mapsto x^k$$ is a bijection. This is equivalent to saying that the equation $x^k = 1$ has a unique solution in $\mathbb{F}_p$. According to this answer here, that is equivalent to $gcd(k,p-1) = 1$. That means that you can take unique square roots only over $\mathbb{F}_2$, over all other finite prime fields there are elements that do not have a square root.
So it turns out that the case of $\mathbb{F}_2$ is the easy one here, if you assume that such a matrix $D$ exists. For very general $M$ it seems finding a root is NP complete, so I don't think you will find an efficient way to do it.