diff --git a/qiskit_optimization/problems/quadratic_expression.py b/qiskit_optimization/problems/quadratic_expression.py
index 2e1501260..c6a7b8a7d 100644
--- a/qiskit_optimization/problems/quadratic_expression.py
+++ b/qiskit_optimization/problems/quadratic_expression.py
@@ -1,6 +1,6 @@
 # This code is part of a Qiskit project.
 #
-# (C) Copyright IBM 2019, 2023.
+# (C) Copyright IBM 2019, 2024.
 #
 # This code is licensed under the Apache License, Version 2.0. You may
 # obtain a copy of this license in the LICENSE.txt file in the root directory
@@ -14,6 +14,7 @@
 
 from typing import List, Union, Dict, Tuple, Any
 
+import sys
 import numpy as np
 from numpy import ndarray
 from scipy.sparse import spmatrix, dok_matrix, tril, triu
@@ -105,6 +106,11 @@ def _coeffs_to_dok_matrix(
         if isinstance(coefficients, (list, ndarray, spmatrix)):
             coefficients = dok_matrix(coefficients)
         elif isinstance(coefficients, dict):
+
+            # Check if the python version is at least 3.9.
+            # See pull request #594 for more details why we do this check.
+            new_update_rule = sys.version_info >= (3, 9)
+
             n = self.quadratic_program.get_num_vars()
             coeffs = dok_matrix((n, n))
             for (i, j), value in coefficients.items():
@@ -112,8 +118,19 @@ def _coeffs_to_dok_matrix(
                     i = self.quadratic_program.variables_index[i]
                 if isinstance(j, str):
                     j = self.quadratic_program.variables_index[j]
-                coeffs[i, j] = value
-            coefficients = coeffs
+
+                if new_update_rule:
+                    if i > j:
+                        coeffs[j, i] += value
+                    else:
+                        coeffs[i, j] += value
+                else:
+                    coeffs[i, j] = value
+
+            if new_update_rule:
+                return coeffs
+            else:
+                coefficients = coeffs
         else:
             raise QiskitOptimizationError(f"Unsupported format for coefficients: {coefficients}")
         return self._triangle_matrix(coefficients)
diff --git a/test/problems/test_quadratic_expression.py b/test/problems/test_quadratic_expression.py
index da3b18c62..394bc94c9 100644
--- a/test/problems/test_quadratic_expression.py
+++ b/test/problems/test_quadratic_expression.py
@@ -1,6 +1,6 @@
 # This code is part of a Qiskit project.
 #
-# (C) Copyright IBM 2020, 2023.
+# (C) Copyright IBM 2020, 2024.
 #
 # This code is licensed under the Apache License, Version 2.0. You may
 # obtain a copy of this license in the LICENSE.txt file in the root directory
@@ -273,6 +273,38 @@ def test_str_repr(self):
             self.assertEqual(str(expr), expected)
             self.assertEqual(repr(expr), f"<QuadraticExpression: {expected[:DEFAULT_TRUNCATE]}...>")
 
+    def test_coeffs_to_dok_matrix(self):
+        """test that the conversion from coefficients to a dok_matrix is correct"""
+        num_vars = 4
+        quadratic_program = QuadraticProgram()
+        for _ in range(num_vars):
+            quadratic_program.binary_var()
+
+        coefficients = {
+            (1, 1): 1,
+            (1, 2): 2,
+            (2, 1): 3,
+            (1, 3): 3,
+            (3, 1): 3,
+            (2, 2): 4,
+            (2, 3): 12,
+            (3, 3): 9,
+        }
+
+        quadratic = QuadraticExpression(quadratic_program, coefficients)
+
+        res = quadratic._coeffs_to_dok_matrix(coefficients)
+
+        expected = dok_matrix((num_vars, num_vars))
+        expected[1, 1] = 1
+        expected[1, 2] = 5
+        expected[1, 3] = 6
+        expected[2, 2] = 4
+        expected[2, 3] = 12
+        expected[3, 3] = 9
+
+        self.assertTrue(np.allclose(expected.todense(), res.todense()))
+
 
 if __name__ == "__main__":
     unittest.main()