diff --git a/CHANGELOG.md b/CHANGELOG.md
index a141d69c7..0d4fe9d67 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -18,6 +18,7 @@
 - Fixed lotsizing_lazy example
 - Fixed incorrect getVal() result when _bestSol.sol was outdated
 - Fixed segmentation fault when using Variable or Constraint objects after freeTransform() or Model destruction
+- getTermsQuadratic() now correctly returns all linear terms
 ### Changed
 - changed default value of enablepricing flag to True
 - Speed up MatrixExpr.sum(axis=...) via quicksum
diff --git a/src/pyscipopt/scip.pxi b/src/pyscipopt/scip.pxi
index b293ed2e2..9c21942bd 100644
--- a/src/pyscipopt/scip.pxi
+++ b/src/pyscipopt/scip.pxi
@@ -8309,8 +8309,16 @@ cdef class Model:
         Returns
         -------
         bilinterms : list of tuple
+            Triples ``(var1, var2, coef)`` for terms of the form
+            ``coef * var1 * var2`` with ``var1 != var2``.
         quadterms : list of tuple
+            Triples ``(var, sqrcoef, lincoef)`` for variables that appear in
+            quadratic or bilinear terms. ``sqrcoef`` is the coefficient of
+            ``var**2``, and ``lincoef`` is the linear coefficient of ``var``
+            if it also appears linearly.
         linterms : list of tuple
+            Pairs ``(var, coef)`` for purely linear variables, i.e.,
+            variables that do not participate in any quadratic or bilinear term.
 
         """
         cdef SCIP_EXPR* expr
@@ -8329,6 +8337,7 @@ cdef class Model:
         cdef int nbilinterms
 
         # quadratic terms
+        cdef SCIP_EXPR* quadexpr
         cdef SCIP_EXPR* sqrexpr
         cdef SCIP_Real sqrcoef
         cdef int nquadterms
@@ -8341,16 +8350,20 @@ cdef class Model:
         assert self.checkQuadraticNonlinear(cons), "constraint is not quadratic"
 
         expr = SCIPgetExprNonlinear(cons.scip_cons)
-        SCIPexprGetQuadraticData(expr, NULL, &nlinvars, &linexprs, &lincoefs, &nquadterms, &nbilinterms, NULL, NULL)
+        SCIPexprGetQuadraticData(expr, NULL, &nlinvars, &linexprs, &lincoefs,
+                                 &nquadterms, &nbilinterms, NULL, NULL)
 
         linterms   = []
         bilinterms = []
-        quadterms  = []
 
+        # Purely linear terms (variables not in any quadratic/bilinear term)
         for termidx in range(nlinvars):
             var = self._getOrCreateVar(SCIPgetVarExprVar(linexprs[termidx]))
             linterms.append((var, lincoefs[termidx]))
 
+        # Collect quadratic terms in a dict so we can merge entries for the same variable.
+        quaddict = {}  # var.ptr() -> [var, sqrcoef, lincoef]
+
         for termidx in range(nbilinterms):
             SCIPexprGetQuadraticBilinTerm(expr, termidx, &bilinterm1, &bilinterm2, &bilincoef, NULL, NULL)
             scipvar1 = SCIPgetVarExprVar(bilinterm1)
@@ -8358,16 +8371,28 @@ cdef class Model:
             var1 = self._getOrCreateVar(scipvar1)
             var2 = self._getOrCreateVar(scipvar2)
             if scipvar1 != scipvar2:
-                bilinterms.append((var1,var2,bilincoef))
+                bilinterms.append((var1, var2, bilincoef))
             else:
-                quadterms.append((var1,bilincoef,0.0))
-
+                # Squared term reported as bilinear var*var
+                key = var1.ptr()
+                if key in quaddict:
+                    quaddict[key][1] += bilincoef
+                else: # TODO: SCIP handles expr like x**2 appropriately, but PySCIPOpt requires this. Need to investigate why.
+                    quaddict[key] = [var1, bilincoef, 0.0]
+
+        # Also collect linear coefficients from the quadratic terms
         for termidx in range(nquadterms):
-            SCIPexprGetQuadraticQuadTerm(expr, termidx, NULL, &lincoef, &sqrcoef, NULL, NULL, &sqrexpr)
-            if sqrexpr == NULL:
-                continue
-            var = self._getOrCreateVar(SCIPgetVarExprVar(sqrexpr))
-            quadterms.append((var,sqrcoef,lincoef))
+            SCIPexprGetQuadraticQuadTerm(expr, termidx, &quadexpr, &lincoef, &sqrcoef, NULL, NULL, &sqrexpr)
+            scipvar1 = SCIPgetVarExprVar(quadexpr)
+            var = self._getOrCreateVar(scipvar1)
+            key = var.ptr()
+            if key in quaddict:
+                quaddict[key][1] += sqrcoef
+                quaddict[key][2] += lincoef
+            else:
+                quaddict[key] = [var, sqrcoef, lincoef]
+
+        quadterms = [tuple(entry) for entry in quaddict.values()]
 
         return (bilinterms, quadterms, linterms)
 
diff --git a/tests/test_nonlinear.py b/tests/test_nonlinear.py
index 383532f2e..55fe6900c 100644
--- a/tests/test_nonlinear.py
+++ b/tests/test_nonlinear.py
@@ -288,6 +288,53 @@ def test_quad_coeffs():
     assert linterms[0][0].name == z.name
     assert linterms[0][1] == 4
 
+
+def test_quad_coeffs_mixed_linear_and_quadratic():
+
+    scip = Model()
+
+    var1 = scip.addVar(name="var1", vtype='C', lb=None)
+    var2 = scip.addVar(name="var2", vtype='C')
+    var3 = scip.addVar(name="var3", vtype='B')
+    var4 = scip.addVar(name="var4", vtype='B')
+
+    cons = scip.addCons(
+        8 * var4
+        + 4 * var3
+        - 5 * var2
+        + 6 * var3 ** 2
+        - 3 * var1 ** 2
+        + 2 * var2 * var1
+        + 7 * var1 * var3
+        == -2
+    )
+
+    bilinterms, quadterms, linterms = scip.getTermsQuadratic(cons)
+
+    # linterms contains only purely linear variables (not in any quadratic/bilinear term)
+    lin_only = {v.name: c for (v, c) in linterms}
+    assert lin_only["var4"] == 8
+    assert len(linterms) == 1  # only var4 is purely linear
+
+    # quadterms contains all variables that appear in quadratic/bilinear terms,
+    # with both their squared coefficient and linear coefficient
+    quad_dict = {v.name: (sqrcoef, lincoef) for v, sqrcoef, lincoef in quadterms}
+    assert quad_dict["var3"] == (6.0, 4.0)   # 6*var3^2 + 4*var3
+    assert quad_dict["var1"] == (-3.0, 0.0)  # -3*var1^2, no linear term
+    assert quad_dict["var2"] == (0.0, -5.0)  # -5*var2, no squared term
+
+    # Verify we can reconstruct all linear coefficients by combining linterms and quadterms
+    full_lin = {}
+    for v, c in linterms:
+        full_lin[v.name] = full_lin.get(v.name, 0.0) + c
+    for v, _, lincoef in quadterms:
+        if lincoef != 0.0:
+            full_lin[v.name] = full_lin.get(v.name, 0.0) + lincoef
+
+    assert full_lin["var4"] == 8
+    assert full_lin["var3"] == 4
+    assert full_lin["var2"] == -5
+
 def test_addExprNonLinear():
     m = Model()
     x = m.addVar("x", lb=0, ub=1, obj=10)