diff --git a/src/code.c b/src/code.c
index 8e8ddb6794f..89bed6a5415 100644
--- a/src/code.c
+++ b/src/code.c
@@ -2450,8 +2450,13 @@ void CodeElmListLevel(CodeState * cs, Int narg, UInt level)
 {
     Expr                ref;            // reference, result
 
+    GAP_ASSERT(narg == 1 || narg == 2);
+
     // allocate the reference and enter the level
-    ref = NewExpr(cs, EXPR_ELM_LIST_LEV, (narg + 2) * sizeof(Expr));
+    if (narg == 1)
+      ref = NewExpr(cs, EXPR_ELM_LIST_LEV, 3 * sizeof(Expr));
+    else // if (narg == 2)
+      ref = NewExpr(cs, EXPR_ELM_MAT_LEV, 4 * sizeof(Expr));
     WRITE_EXPR(cs, ref, narg + 1, level);
 
     // let 'CodeElmListUniv' do the rest
diff --git a/src/code.h b/src/code.h
index a66c7f40a63..3bd16790381 100644
--- a/src/code.h
+++ b/src/code.h
@@ -534,6 +534,7 @@ enum EXPR_TNUM {
     EXPR_ELM_MAT,
     EXPR_ELMS_LIST,
     EXPR_ELM_LIST_LEV,
+    EXPR_ELM_MAT_LEV,
     EXPR_ELMS_LIST_LEV,
     EXPR_ISB_LIST,
 
diff --git a/src/syntaxtree.c b/src/syntaxtree.c
index 8184c60a25f..5838318efed 100644
--- a/src/syntaxtree.c
+++ b/src/syntaxtree.c
@@ -904,9 +904,14 @@ static const CompilerT Compilers[] = {
     COMPILER_(
         EXPR_ELM_MAT, ARG_EXPR_("list"), ARG_EXPR_("row"), ARG_EXPR_("col")),
     COMPILER_(EXPR_ELMS_LIST, ARG_EXPR_("list"), ARG_EXPR_("poss")),
+    COMPILER_(EXPR_ELM_MAT_LEV,
+              ARG_EXPR_("matrices"),
+              ARG_EXPR_("row"),
+              ARG_EXPR_("col"),
+              ARG_EXPR("level", ObjInt_UInt, SyntaxTreeCodeObjInt)),
     COMPILER_(EXPR_ELM_LIST_LEV,
               ARG_EXPR_("lists"),
-              ARGS_EXPR("pos"),
+              ARG_EXPR_("pos"),
               ARG_EXPR("level", ObjInt_UInt, SyntaxTreeCodeObjInt)),
     COMPILER_(EXPR_ELMS_LIST_LEV,
               ARG_EXPR_("lists"),
diff --git a/src/vars.c b/src/vars.c
index c15583e032b..95f2471fd61 100644
--- a/src/vars.c
+++ b/src/vars.c
@@ -825,24 +825,19 @@ static Obj EvalElmListLevel(Expr expr)
 {
     Obj                 lists;          // lists, left operand
     Obj                 pos;            // position, right operand
-    Obj                 ixs;
     UInt                level;          // level
-    Int narg;
-    Int i;
+    Obj                 ixs;
 
     // evaluate lists (if this works, then <lists> is nested <level> deep,
     // checking it is nested <level>+1 deep is done by 'ElmListLevel')
     lists = EVAL_EXPR(READ_EXPR(expr, 0));
-    narg = SIZE_EXPR(expr)/sizeof(Expr) -2;
-    ixs = NEW_PLIST(T_PLIST, narg);
-    for (i = 1; i <= narg; i++) {
-      pos = EVAL_EXPR( READ_EXPR(expr, i));
-      SET_ELM_PLIST(ixs, i, pos);
-      CHANGED_BAG(ixs);
-    }
-    SET_LEN_PLIST(ixs, narg);
-    // get the level
-    level = READ_EXPR(expr, narg + 1);
+    pos = EVAL_EXPR(READ_EXPR(expr, 1));
+    level = READ_EXPR(expr, 2);
+
+    ixs = NEW_PLIST(T_PLIST, 1);
+    SET_ELM_PLIST(ixs, 1, pos);
+    CHANGED_BAG(ixs);
+    SET_LEN_PLIST(ixs, 1);
 
     // select the elements from several lists (store them in <lists>)
     ElmListLevel( lists, ixs, level );
@@ -852,6 +847,49 @@ static Obj EvalElmListLevel(Expr expr)
 }
 
 
+/****************************************************************************
+**
+*F  EvalElmMatLevel(<expr>)  . . . . . .  select elements of several matrices
+**
+**  'EvalElmMatLevel' evaluates the matrix element  expression <expr> of  the
+**  form '<list>...{<positions>}...[<row>,<col>]',  where there  may actually
+**  be several '{<positions>}' selections between <list> and  '[<position>]'.
+**  The  number of those is called   the level.  'EvalElmMatLevel'  goes that
+**  deep  into the left operand  and  selects the element at <row>,<col> from
+**  each of those matrices. For example,  if the level is 1, the left operand
+**  must be a list of matrices and 'EvalElmMatLevel'  selects the element  at
+**  <row>,<col>  from each  of  the matrices  and returns  the list of  those
+**  values.
+*/
+static Obj EvalElmMatLevel(Expr expr)
+{
+    Obj                 matrices;       // matrices
+    Obj                 row;            // row position
+    Obj                 col;            // column position
+    UInt                level;          // level
+    Obj                 ixs;
+
+    // evaluate matrices (if this works, then <matrices> is nested <level>
+    // deep, checking it is nested <level>+1 deep is done by 'ElmListLevel')
+    matrices = EVAL_EXPR(READ_EXPR(expr, 0));
+    row = EVAL_EXPR(READ_EXPR(expr, 1));
+    col = EVAL_EXPR(READ_EXPR(expr, 2));
+    level = READ_EXPR(expr, 3);
+
+    ixs = NEW_PLIST(T_PLIST, 2);
+    SET_ELM_PLIST(ixs, 1, row);
+    SET_ELM_PLIST(ixs, 2, col);
+    CHANGED_BAG(ixs);
+    SET_LEN_PLIST(ixs, 2);
+
+    // select the elements from several matrices (store them in <matrices>)
+    ElmListLevel( matrices, ixs, level );
+
+    // return the elements
+    return matrices;
+}
+
+
 /****************************************************************************
 **
 *F  EvalElmsListLevel(<expr>) . . .  select several elements of several lists
@@ -1036,16 +1074,21 @@ static void PrintElmMat(Expr expr)
 
 static void PrintElmListLevel(Expr expr)
 {
-    Int i;
-    Int narg = SIZE_EXPR(expr)/sizeof(Expr) -2 ;
     Pr("%2>", 0, 0);
     PrintExpr(READ_EXPR(expr, 0));
     Pr("%<[", 0, 0);
     PrintExpr(READ_EXPR(expr, 1));
-    for (i = 2; i <= narg; i++) {
-        Pr("%<, %>", 0, 0);
-        PrintExpr(READ_EXPR(expr, i));
-    }
+    Pr("%<]", 0, 0);
+}
+
+static void PrintElmMatLevel(Expr expr)
+{
+    Pr("%2>", 0, 0);
+    PrintExpr(READ_EXPR(expr, 0));
+    Pr("%<[", 0, 0);
+    PrintExpr(READ_EXPR(expr, 1));
+    Pr("%<, %>", 0, 0);
+    PrintExpr(READ_EXPR(expr, 2));
     Pr("%<]", 0, 0);
 }
 
@@ -2248,8 +2291,10 @@ static Int InitKernel (
     // install executors, evaluators, and printers for matrix elements
     InstallExecStatFunc(STAT_ASS_MAT, ExecAssMat);
     InstallEvalExprFunc(EXPR_ELM_MAT, EvalElmMat);
+    InstallEvalExprFunc(EXPR_ELM_MAT_LEV, EvalElmMatLevel);
     InstallPrintStatFunc(STAT_ASS_MAT, PrintAssMat);
     InstallPrintExprFunc(EXPR_ELM_MAT, PrintElmMat);
+    InstallPrintExprFunc(EXPR_ELM_MAT_LEV, PrintElmMatLevel);
 
     // install executors, evaluators, and printers for record elements
     InstallExecStatFunc( STAT_ASS_REC_NAME   , ExecAssRecName);