diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8a4866b00..aeab43ace 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,7 @@
 - Added stage checks to presolve, freereoptsolve, freetransform
 - Added primal_dual_evolution recipe and a plot recipe
 ### Fixed
+- Only redirect stdout and stderr in redirectOutput() so that file output still works afterwards
 ### Changed
 - GitHub actions using Mac now use precompiled SCIP from latest release
 ### Removed
diff --git a/src/pyscipopt/scip.pxi b/src/pyscipopt/scip.pxi
index 04dba32f2..bf575ba65 100644
--- a/src/pyscipopt/scip.pxi
+++ b/src/pyscipopt/scip.pxi
@@ -12,7 +12,7 @@ cimport cython
 from cpython cimport Py_INCREF, Py_DECREF
 from cpython.pycapsule cimport PyCapsule_New, PyCapsule_IsValid, PyCapsule_GetPointer
 from libc.stdlib cimport malloc, free
-from libc.stdio cimport fdopen, fclose
+from libc.stdio cimport stdout, stderr, fdopen, fputs, fflush, fclose
 from posix.stdio cimport fileno
 
 from collections.abc import Iterable
@@ -1851,10 +1851,22 @@ cdef class Constraint:
 
 
 cdef void relayMessage(SCIP_MESSAGEHDLR *messagehdlr, FILE *file, const char *msg) noexcept:
-    sys.stdout.write(msg.decode('UTF-8'))
+    if file is stdout:
+        sys.stdout.write(msg.decode('UTF-8'))
+    elif file is stderr:
+        sys.stderr.write(msg.decode('UTF-8'))
+    else:
+        if msg is not NULL:
+            fputs(msg, file)
+        fflush(file)
 
 cdef void relayErrorMessage(void *messagehdlr, FILE *file, const char *msg) noexcept:
-    sys.stderr.write(msg.decode('UTF-8'))
+    if file is NULL:
+        sys.stderr.write(msg.decode('UTF-8'))
+    else:
+        if msg is not NULL:
+            fputs(msg, file)
+        fflush(file)
 
 # - remove create(), includeDefaultPlugins(), createProbBasic() methods
 # - replace free() by "destructor"
diff --git a/tests/test_model.py b/tests/test_model.py
index 2d07331e5..a87548e99 100644
--- a/tests/test_model.py
+++ b/tests/test_model.py
@@ -3,6 +3,7 @@
 import itertools
 
 from pyscipopt import Model, SCIP_STAGE, SCIP_PARAMSETTING, quicksum
+from helpers.utils import random_mip_1
 
 def test_model():
     # create solver instance
@@ -477,4 +478,36 @@ def test_getObjVal():
     assert m.getVal(x) == 0
 
     assert m.getObjVal() == 16 
-    assert m.getVal(x) == 0
\ No newline at end of file
+    assert m.getVal(x) == 0
+
+# tests writeProblem() after redirectOutput()
+def test_redirection():
+
+    # create problem instances
+    original = random_mip_1(False, False, False, -1, True)
+    redirect = Model()
+
+    # redirect console output
+    original.redirectOutput()
+
+    # write problem instance
+    original.writeProblem("redirection.lp")
+
+    # solve original instance
+    original.optimize()
+
+    # read problem instance
+    redirect.readProblem("redirection.lp")
+
+    # remove problem file
+    os.remove("redirection.lp")
+
+    # compare problem dimensions
+    assert redirect.getNVars(False) == original.getNVars(False)
+    assert redirect.getNConss(False) == original.getNConss(False)
+
+    # solve redirect instance
+    redirect.optimize()
+
+    # compare objective values
+    assert original.isEQ(redirect.getObjVal(), original.getObjVal())