diff --git a/CHANGELOG.md b/CHANGELOG.md
index addf21992..72010d17a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -18,6 +18,7 @@
 ### Fixed
 - Raised an error when an expression is used when a variable is required
 - Fixed some compile warnings
+- _VarArray can accept MatrixVariable now
 ### Changed
 - MatrixExpr.sum() now supports axis arguments and can return either a scalar or MatrixExpr, depending on the result dimensions.
 - AddMatrixCons() also accepts ExprCons.
diff --git a/src/pyscipopt/scip.pxi b/src/pyscipopt/scip.pxi
index a380b1456..66d14426b 100644
--- a/src/pyscipopt/scip.pxi
+++ b/src/pyscipopt/scip.pxi
@@ -2499,21 +2499,18 @@ cdef class _VarArray:
 
     def __cinit__(self, object vars):
         if isinstance(vars, Variable):
-            self.size = 1
-            self.ptr = <SCIP_VAR**> malloc(sizeof(SCIP_VAR*))
-            self.ptr[0] = (<Variable>vars).scip_var
+            vars = [vars]
+        elif isinstance(vars, (list, tuple, MatrixVariable)):
+            vars = np.ravel(vars)
         else:
-            if not isinstance(vars, (list, tuple)):
-                raise TypeError("Expected Variable or list of Variable, got %s." % type(vars))
-            self.size = len(vars)
-            if self.size == 0:
-                self.ptr = NULL
-            else:
-                self.ptr = <SCIP_VAR**> malloc(self.size * sizeof(SCIP_VAR*))
-                for i, var in enumerate(vars):
-                    if not isinstance(var, Variable):
-                        raise TypeError("Expected Variable, got %s." % type(var))
-                    self.ptr[i] = (<Variable>var).scip_var
+            raise TypeError(f"Expected Variable or list of Variable, got {type(vars)}.")
+
+        self.size = len(vars)
+        self.ptr = <SCIP_VAR**> malloc(self.size * sizeof(SCIP_VAR*)) if self.size else NULL
+        for i, var in enumerate(vars):
+            if not isinstance(var, Variable):
+                raise TypeError(f"Expected Variable, got {type(var)}.")
+            self.ptr[i] = (<Variable>var).scip_var
 
     def __dealloc__(self):
         if self.ptr != NULL:
diff --git a/tests/test_cons.py b/tests/test_cons.py
index 36a0dfcec..756088612 100644
--- a/tests/test_cons.py
+++ b/tests/test_cons.py
@@ -179,6 +179,40 @@ def test_cons_indicator():
     assert m.isEQ(m.getVal(x), 1)
     assert c1.getConshdlrName() == "indicator"
 
+def test_cons_indicator_with_matrix_binvar():
+    # test matrix variable binvar #1043
+    m = Model()
+
+    with pytest.raises(TypeError):
+        m.addConsIndicator(m.addVar(vtype="B") <= 1, 1)
+
+    # test binvar with (1, 1, 1) shape of matrix variable
+    x = m.addVar(vtype="B")
+    binvar1 = m.addMatrixVar(((1, 1, 1)), vtype="B")
+    m.addConsIndicator(x >= 1, binvar1, activeone=True)
+    m.addConsIndicator(x <= 0, binvar1, activeone=False)
+
+    # test binvar with (2, 3) shape of matrix variable
+    y = m.addVar(vtype="B")
+    binvar2 = m.addMatrixVar(((2, 3)), vtype="B")
+    m.addConsIndicator(y >= 1, binvar2, activeone=True)
+    m.addConsIndicator(y <= 0, binvar2, activeone=False)
+
+    # test binvar with (2, 1) shape of list of lists
+    z = m.addVar(vtype="B")
+    binvar3 = [[m.addVar(vtype="B")], [m.addVar(vtype="B")]]
+    m.addConsIndicator(z >= 1, binvar3, activeone=True)
+    m.addConsIndicator(z <= 0, binvar3, activeone=False)
+
+    m.setObjective(
+        binvar1.sum() + binvar2.sum() + binvar3[0][0] + binvar3[1][0], "maximize"
+    )
+    m.optimize()
+
+    assert m.getVal(x) == 1
+    assert m.getVal(y) == 1
+    assert m.getVal(z) == 1
+
 @pytest.mark.xfail(
     reason="addConsIndicator doesn't behave as expected when binary variable is False. See Issue #717."
 )