diff --git a/CHANGELOG.md b/CHANGELOG.md
index c19ccfa4f..70b43499a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Unreleased
 ### Added
+- Added `getCountedSols()` to retrieve the solutions collected by `count()` (#906)
 ### Fixed
 ### Changed
 ### Removed
diff --git a/src/pyscipopt/scip.pxd b/src/pyscipopt/scip.pxd
index fe5f3ac6b..bf5479552 100644
--- a/src/pyscipopt/scip.pxd
+++ b/src/pyscipopt/scip.pxd
@@ -522,6 +522,9 @@ cdef extern from "scip/scip.h":
     ctypedef struct SCIP_DIVESET:
         pass
 
+    ctypedef struct SCIP_SPARSESOL:
+        pass
+
     ctypedef struct SCIP_HASHMAP:
         pass
 
@@ -2102,6 +2105,9 @@ cdef extern from "scip/cons_countsols.h":
     SCIP_RETCODE SCIPcount(SCIP* scip)
     SCIP_RETCODE SCIPsetParamsCountsols(SCIP* scip)
     SCIP_Longint SCIPgetNCountedSols(SCIP* scip, SCIP_Bool* valid)
+    void SCIPgetCountedSparseSols(SCIP* scip, SCIP_VAR*** vars, int* nvars, SCIP_SPARSESOL*** sols, int* nsols)
+    void SCIPsparseSolGetFirstSol(SCIP_SPARSESOL* sparsesol, SCIP_Longint* sol, int nvars)
+    SCIP_Bool SCIPsparseSolGetNextSol(SCIP_SPARSESOL* sparsesol, SCIP_Longint* sol, int nvars)
 
 cdef extern from "scip/paramset.h":
 
diff --git a/src/pyscipopt/scip.pxi b/src/pyscipopt/scip.pxi
index 52123c2aa..f7b449022 100644
--- a/src/pyscipopt/scip.pxi
+++ b/src/pyscipopt/scip.pxi
@@ -12545,6 +12545,55 @@ cdef class Model:
             print('total number of solutions found is not valid!')
         return nsols
 
+    def getCountedSols(self):
+        """
+        Get the feasible solutions collected by `count()`.
+
+        Requires `constraints/countsols/collect` to be True and presolving to be off
+        before calling `count()`. Only integral variables are counted, so continuous
+        variables are not part of the solutions. The solutions are freed together with
+        the transformed problem.
+
+        Returns
+        -------
+        list of dict
+            Each solution maps the names of the transformed variables to their values.
+
+        """
+        cdef SCIP_VAR** vars
+        cdef SCIP_SPARSESOL** sparsesols
+        cdef SCIP_Longint* sol
+        cdef SCIP_Bool valid
+        cdef int nvars
+        cdef int nsparsesols
+        cdef int i
+        cdef int s
+
+        if not self.getParam("constraints/countsols/collect"):
+            raise Warning("set constraints/countsols/collect to True before calling count()")
+        if self.getParam("presolving/maxrounds") != 0:
+            raise Warning("turn presolving off before calling count()")
+
+        SCIPgetNCountedSols(self._scip, &valid)
+        if not valid:
+            raise Warning("number of counted solutions does not fit into a long integer, cannot collect them")
+
+        SCIPgetCountedSparseSols(self._scip, &vars, &nvars, &sparsesols, &nsparsesols)
+        names = [bytes(SCIPvarGetName(vars[i])).decode("utf-8") for i in range(nvars)]
+
+        sol = <SCIP_Longint*> malloc(nvars * sizeof(SCIP_Longint))
+        sols = []
+        for s in range(nsparsesols):
+            SCIPsparseSolGetFirstSol(sparsesols[s], sol, nvars)
+            while True:
+                sols.append({names[i]: sol[i] for i in range(nvars)})
+                if not SCIPsparseSolGetNextSol(sparsesols[s], sol, nvars):
+                    break
+
+        free(sol)
+
+        return sols
+
     def setParamsCountsols(self):
         """Sets SCIP parameters such that a valid counting process is possible."""
         PY_SCIP_CALL(SCIPsetParamsCountsols(self._scip))
diff --git a/src/pyscipopt/scip.pyi b/src/pyscipopt/scip.pyi
index 316b13e9f..1cc29b1a4 100644
--- a/src/pyscipopt/scip.pyi
+++ b/src/pyscipopt/scip.pyi
@@ -1226,6 +1226,7 @@ class Model:
     def getConsVars(self, constraint: Constraint) -> list[Variable] | None: ...
     def getConss(self, transformed: bool = True) -> list[Constraint]: ...
     def getFocusNode(self) -> Node | None: ...
+    def getCountedSols(self) -> list[dict[str, int]]: ...
     def getCurrentNode(self) -> Node | None: ...
     def getCutEfficacy(self, cut: Row, sol: Solution | None = None) -> float: ...
     def getCutLPSolCutoffDistance(self, cut: Row, sol: Solution) -> float: ...
diff --git a/tests/test_model.py b/tests/test_model.py
index 6736fbb5f..c4ee90b0e 100644
--- a/tests/test_model.py
+++ b/tests/test_model.py
@@ -662,3 +662,40 @@ def test_memory_methods():
 
     assert m.getMemUsed() > 0
     assert m.getMemTotal() > 0
+
+
+def test_getCountedSols():
+    m = Model()
+    m.hideOutput()
+    m.setParam("constraints/countsols/collect", True)
+    m.setPresolve(SCIP_PARAMSETTING.OFF)
+    x = m.addVar("x", vtype="I", lb=-2, ub=3)
+    y = m.addVar("y", vtype="I", lb=0, ub=4)
+    z = m.addVar("z", vtype="I", lb=0, ub=10)
+    w = m.addVar("w", vtype="B")
+    m.addVar("c", lb=0, ub=1)
+    m.addCons(z == 2*x + 3)
+    m.addCons(x + y <= 4)
+    m.addCons(w + y >= 1)
+
+    assert m.getCountedSols() == []
+
+    m.count()
+    sols = m.getCountedSols()
+    expected = {(a, b, c, d) for a, b, c, d in itertools.product(range(-2, 4), range(5), range(11), range(2))
+                if c == 2*a + 3 and a + b <= 4 and d + b >= 1}
+    assert len(sols) == m.getNCountedSols() == len(expected)
+    assert all(set(sol) == {"t_x", "t_y", "t_z", "t_w"} for sol in sols)
+    assert {tuple(sol["t_" + k] for k in "xyzw") for sol in sols} == expected
+
+    m.freeTransform()
+    assert m.getCountedSols() == []
+
+    m.setParam("constraints/countsols/collect", False)
+    with pytest.raises(Warning):
+        m.getCountedSols()
+
+    m.setParam("constraints/countsols/collect", True)
+    m.setPresolve(SCIP_PARAMSETTING.DEFAULT)
+    with pytest.raises(Warning):
+        m.getCountedSols()