In order to answer the question we need a finite presentation
of ${\rm SL}(3,\mathbb{Z})$ and a general method to find all subgroups of index
$\leq n$ of a finitely presented group:
A finite presentation for ${\rm SL}(3,\mathbb{Z})$ can be found for example in Theorem 2 in
Marston Conder, Edmund Robertson, Peter Williams:
Presentations for 3-dimensional special linear groups over integer rings,
Proc. Amer. Math. Soc. 115 (1992), no. 1, 19-26.
http://www.ams.org/journals/proc/1992-115-01/S0002-9939-1992-1079696-5/S0002-9939-1992-1079696-5.pdf.
The finite presentation given in this paper is
$$ {\rm SL}(3,\mathbb{Z}) \cong \left< x, y, z \ | \ x^3 = y^3 = z^2 = (xz)^3 = (yz)^3 = (x^{-1}zxy)^2 = (y^{-1}zyx)^2 = (xy)^6 = 1 \right>
$$
on the generators
$$
x \ = \
\left(
\begin{array}{rrr}
1 & 0 & 1 \\\
0 & -1 & -1 \\\
0 & 1 & 0
\end{array}
\right), \ \
y \ = \
\left(
\begin{array}{rrr}
0 & 1 & 0 \\\
0 & 0 & 1 \\\
1 & 0 & 0
\end{array}
\right), \ \
z \ = \
\left(
\begin{array}{rrr}
0 & 1 & 0 \\\
1 & 0 & 0 \\\
-1 & -1 & -1
\end{array}
\right).
$$
A general method to find all subgroups of index $\leq n$ of a finitely presented
group is the so-called low index subgroups procedure.
This algorithm is described in Section 5.4 in
Derek F. Holt, Bettina Eick, and Eamonn A. O'Brien, Handbook of
computational group theory, Discrete Mathematics and its Applications (Boca
Raton), Chapman & Hall / CRC, Boca Raton, FL, 2005. MR 2129747 (2006f:20001).
For an online resource, see e.g.
Marston Conder: Applications and adaptations of the low index subgroups
procedure, http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.107.5164.
The low index subgroups procedure is implemented in the GAP computer algebra system
(cf. http://www.gap-system.org). Hence all we need to do is to enter the presentation
of ${\rm SL}(3,\mathbb{Z})$ taken from the above paper into GAP ...
gap> F := FreeGroup("x","y","z");;
gap> AssignGeneratorVariables(F);
#I Assigned the global variables [ x, y, z ]
gap> G := F/[x^3,y^3,z^2,(x*z)^3,(y*z)^3,(x^-1*z*x*y)^2,(y^-1*z*y*x)^2,(x*y)^6];
<fp group on the generators [ x, y, z ]>
... and to run the algorithm on it:
gap> sub := LowIndexSubgroupsFpGroup(G,7);;
gap> List(sub,H->Index(G,H));
[ 1, 7, 7 ]
gap> gens := List(sub,GeneratorsOfGroup);
[ [ x, y, z ], [ x, z, y*z*y^-1, (y*x)^2*y ],
[ x, y*x^-1*z^-1, y^-1*z*y, z*y^-1*x*y, y^-1*x*y*x^-1*y ] ]
This tells us that the smallest index of a proper subgroup of ${\rm SL}(3,\mathbb{Z})$
is 7, and that there are 2 conjugacy classes of subgroups of index 7.
Now it is straightforward to obtain generators for our subgroups in terms of matrices:
gap> x := [ [ 1, 0, 1 ], [ 0, -1, -1 ], [ 0, 1, 0 ] ];;
gap> y := [ [ 0, 1, 0 ], [ 0, 0, 1 ], [ 1, 0, 0 ] ];;
gap> z := [ [ 0, 1, 0 ], [ 1, 0, 0 ], [ -1, -1, -1 ] ];;
gap> List(gens[2],g->MappedWord(g,GeneratorsOfGroup(G),[x,y,z]));
[ [ [ 1, 0, 1 ], [ 0, -1, -1 ], [ 0, 1, 0 ] ],
[ [ 0, 1, 0 ], [ 1, 0, 0 ], [ -1, -1, -1 ] ],
[ [ 0, 0, 1 ], [ -1, -1, -1 ], [ 1, 0, 0 ] ],
[ [ -1, -1, -1 ], [ 0, 0, 1 ], [ 0, 1, -1 ] ] ]
gap> List(gens[3],g->MappedWord(g,GeneratorsOfGroup(G),[x,y,z]));
[ [ [ 1, 0, 1 ], [ 0, -1, -1 ], [ 0, 1, 0 ] ],
[ [ -1, -1, -1 ], [ 0, 1, 1 ], [ 0, 0, -1 ] ],
[ [ -1, -1, -1 ], [ 0, 0, 1 ], [ 0, 1, 0 ] ],
[ [ 1, 1, 0 ], [ 0, 0, 1 ], [ 0, -1, 0 ] ],
[ [ -1, 0, -1 ], [ 2, 1, 1 ], [ 0, -1, 0 ] ] ]
So our representatives for the conjugacy classes of subgroups
of ${\rm SL}(3,\mathbb{Z})$ of index 7 are
$$
G_{7,1} \ = \
\left<
\left(\begin{array}{rrr}
1&0&1\\\
0&-1&-1\\\
0&1&0
\end{array}\right), \
\left(\begin{array}{rrr}%
0&1&0\\\
1&0&0\\\
-1&-1&-1
\end{array}\right), \\
\left(\begin{array}{rrr}%
0&0&1\\\
-1&-1&-1\\\
1&0&0
\end{array}\right), \
\left(\begin{array}{rrr}%
-1&-1&-1\\\
0&0&1\\\
0&1&-1
\end{array}\right)
\right>
$$
and
$$
G_{7,2} \ = \
\left<
\left(\begin{array}{rrr}%
1&0&1\\\
0&-1&-1\\\
0&1&0
\end{array}\right), \
\left(\begin{array}{rrr}%
-1&-1&-1\\\
0&1&1\\\
0&0&-1
\end{array}\right), \
\left(\begin{array}{rrr}%
-1&-1&-1\\\
0&0&1\\\
0&1&0
\end{array}\right), \\
\left(\begin{array}{rrr}%
1&1&0\\\
0&0&1\\\
0&-1&0
\end{array}\right), \
\left(\begin{array}{rrr}%
-1&0&-1\\\
2&1&1\\\
0&-1&0
\end{array}\right)
\right>.
$$
The computations above take just a few milliseconds. If one is willing to put in a minute or so,
then one can go a bit further and compute representatives for the conjugacy classes of
subgroups of ${\rm SL}(3,\mathbb{Z})$ of index $\leq 30$:
gap> sub := LowIndexSubgroupsFpGroup(G,30);;
gap> List(sub,H->Index(G,H));
[ 1, 8, 7, 28, 14, 13, 7, 13, 28, 26, 14, 26, 28, 28, 28, 24, 21 ]
So we have subgroups of indices 7, 8, 13, 14, 21, 24, 26 and 28, and there are
no proper subgroups of other indices $\leq 30$.
Generators of the subgroups in terms of our generators $x, y, z$ of
${\rm SL}(3,\mathbb{Z})$ can be determined easily as well:
gap> List(sub,GeneratorsOfGroup);
[ [ x, y, z ], [ x, y ], [ x, z, y*z*y^-1, (y*x)^2*y ],
[ x, z, y*z*y^-1, y*x*(y*x^-1)^2*y, (y*x)^2*(y^-1*x)^2*y^-1 ],
[ x, z, (y*x)^2*y ], [ x, y*x^-1*z^-1, y^-1*z*y, z*y^-1*x^-1*y ],
[ x, y*x^-1*z^-1, y^-1*z*y, z*y^-1*x*y, y^-1*x*y*x^-1*y ],
[ x, y*x^-1*z^-1, y^-1*z*y, z*y^-1*x*y, y^-1*(x*y)^2*x^-1*y^-1*x^-1*y ],
[ x, y*x^-1*z^-1, y^-1*z*y, z*y^-1*x*y,
y^-1*(x*y)^2*x^-1*y^-1*x*y*x*y^-1*x^-1*y,
y^-1*(x*y)^2*x^-1*y*x*y^-1*x*y*x^-1*y^-1*x^-1*y ],
[ x, y*x^-1*z^-1, z*y^-1*x^-1*y ],
[ x, y*z*x^-1*y^-1, z*y*x*z^-1, (y*x)^2*y, y^-1*x*y*x^-1*y ],
[ x, y*z*x^-1*y^-1, z*y*x*z^-1, z*y^-1*x*y, y^-1*(x*y)^2*x^-1*y^-1*x^-1*y ],
[ x, y^-1*z*y, (y*x)^2*y, y^-1*x*y*x^-1*y ],
[ x, y^-1*z*y, (y*x)^2*y, z*y^-1*(x^-1*y)^2 ],
[ x, y^-1*z*y, y*x^-1*z*y^-1*x^-1*y^-1, y^-1*x*y*z^-1*x^-1*y^-1 ],
[ y*x^-1, y^-1*x ], [ z, x*z*x^-1, y*z*y^-1, (y*x)^3 ] ]
$\mathrm{SL}(2,\mathbb{Z})$
, or its quotient by scalars the modular group, the congruence subgroup property fails and the subgroup structure is much richer than in higher ranks. But the modular group has a lot of older literature relative to its action on the complex upper half plane, etc. $\endgroup$