Example of matrix $M\in GL_3(\mathbb{Z}/7\mathbb{Z})$ such that $\langle M; +, \cdot \rangle \simeq GF(7^3)$ and the multiplicative order of $M$ is 3 
I would want to make an example of a matrix $M \in GL_3(\mathbb{Z}_7)$ such that $\langle M; +, \cdot \rangle \simeq GF(7^3)$ and the multiplicative order of $M$ is $3$.

Any hints how to do that with a computer algebra system would be appreciated.
 A: While you figure out what the question is supposed to be, here is some advice on using GAP to find an answer:
The naive way of looking for such a matrix is just to walk through the elements of GL(3,7):
First( GL(3,7), m -> Order(m) = 3 and Size(Algebra(GF(7),[m])) = 7^3 );

However, this runs into a performance problem as GL(3,7) is quite large:
gap> Size(GL(3,7));
33784128

Instead, we use the fact that GL(3,7) has relatively few and relatively well understood conjugacy classes:
gap> ccreps := List( ConjugacyClasses( GL( 3, 7 ) ), Representative );;
gap> Size( ccreps );
336

Now we can examine those 336 conjugacy class representatives to see if any satisfy your question:
gap> ms := Filtered( ccreps, m -> Order(m) = 3 and
> Size(Algebra(GF(7),[m])) = 7^3 );; Size( ms );
1
gap> Perform( ms, Display );
 4 . .
 . 1 .
 . . 2

So we see there is only one matrix of order 3 in GL(3,7), up to conjugacy, such that $\langle M,+,\cdot\rangle$ has $7^3$ elements.  Unfortunately, it is isomorphic to $GF(7)^3$ not $GF(7^3)$.
It is quite easy to find an M that generates a copy of $GF(7^3)$: one simply looks for a matrix whose order divides $7^3-1$ but not $7-1$.  In other words, we just want a matrix of order 57:
gap> ms := Filtered( ccreps, m -> Order(m) = 57 );;
gap> Size( ms );
12

We get 12 conjugacy classes, but in fact these form a single rational class:
gap> ForAll( ms, m -> ForAny(PrimeResidues(57),
> k -> IsConjugate(GL(3,7),m,m^k)));
true

Now we want to check the algebras just to make sure:
gap> as:=List( ms, m -> Algebra(GF(7),[m]));;
gap> ForAll( as, IsSimple );
true
gap> ForAll( as, IsCommutative );
true
gap> ForAll( as, a -> Size(a) = 7^3 );
true

A: Actually you can do the generalized version (described in pol's comment) with Mathematica, too. The task at hand is to find a 3x3 matrix $M$ with entries in $GF(p)=\mathbf{Z}/p\mathbf{Z}$ such that $M$ has multiplicative order $q$, and that the algebra it generates is isomorphic to $GF(p^3)$. The multiplicative group of $GF(p^3)$ is cyclic of order $p^3-1$, so for such an $M$ to exist, it is necessary that $q\mid p^3-1$. For the algebra $GF(p)[M]$ to be a field, it is necessary that the minimal polynomial of $M$ is cubic, and that it does not split into lower degree factors over $GF(p)$. The minimal polynomial of a matrix is always a factor of its characteristic polynomial $\chi_M(x)$, so we only need $\chi_M(x)$ to be irreducible over $GF(p)$. Then $\langle M,+,\cdot\rangle=GF(p)[x]\langle \chi_M(x)\rangle$, which is isomorphic to $GF(p^3)$. 
But the constraint on the order of $M$ forces the minimal polynomial to be a factor of $x^q-1$. So we need $x^q-1$ to have an irreducible factor of degree three. Note that when $q\mid p-1$, then the polynomial $x^q-1$ splits into a product of linear factors over $GF(p)$, so we have to avoid that.
Following Jack's suggestion ($57=3\cdot19$) let's try with $p=7$, $q=19$ with Mathematica
In[1]:=Factor[x^19-1,Modulus->7]
Out[1]=(6+x)(6+2x+x^3)(6+3x+3x^2+x^3)(6+x+4x^2+x^3)(6+4x+4x^2+x^3)(6+5x^2+x^3)(6+3x+6x^2+x^3)

So we have no less than six cubic factors (it was actually possible to tell this in advance, but let's skip that). Let's use the first: $6+2x+x^3$. We use a standard recipe of companion matrix to produce a matrix with a given minimal polynomial, and pick
$$
M=\pmatrix{0&0&-6\cr 1&0&-2\cr 0&1&0\cr}.
$$
It is straightforward to verify that $M^{19}\equiv I_3 \pmod 7$. 
As a final check we see that $q=3, p=7$ does not work:
In[2]:=Factor[x^3-1,Modulus->7]
Out[2]=(3+x)(5+x)(6+x)

confirming the findings from my comments and Jack's answer that the third roots of unity in a field of characteristic seven are $4=-3$, $5=-2$, and $1=-6$.
