diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6e203885c..e58da1c71 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,7 @@
 - all fundamental callbacks now raise an error if not implemented
 - Fixed the type of MatrixExpr.sum(axis=...) result from MatrixVariable to MatrixExpr.
 - Updated IIS result in PyiisfinderExec()
+- Model.getVal now supports GenExpr type
 - Fixed lotsizing_lazy example
 - Fixed incorrect getVal() result when _bestSol.sol was outdated
 ### Changed
diff --git a/src/pyscipopt/expr.pxi b/src/pyscipopt/expr.pxi
index e0d35e8b0..e37d14867 100644
--- a/src/pyscipopt/expr.pxi
+++ b/src/pyscipopt/expr.pxi
@@ -42,9 +42,20 @@
 # which should, in princple, modify the expr. However, since we do not implement __isub__, __sub__
 # gets called (I guess) and so a copy is returned.
 # Modifying the expression directly would be a bug, given that the expression might be re-used by the user. </pre>
+import math
+from typing import TYPE_CHECKING
+
+from pyscipopt.scip cimport Variable, Solution
+from cpython.dict cimport PyDict_Next
+from cpython.ref cimport PyObject
+
 import numpy as np
 
 
+if TYPE_CHECKING:
+    double = float
+
+
 def _is_number(e):
     try:
         f = float(e)
@@ -87,23 +98,25 @@ def _expr_richcmp(self, other, op):
         raise NotImplementedError("Can only support constraints with '<=', '>=', or '=='.")
 
 
-class Term:
+cdef class Term:
     '''This is a monomial term'''
 
-    __slots__ = ('vartuple', 'ptrtuple', 'hashval')
+    cdef readonly tuple vartuple
+    cdef readonly tuple ptrtuple
+    cdef Py_ssize_t hashval
 
-    def __init__(self, *vartuple):
+    def __init__(self, *vartuple: Variable):
         self.vartuple = tuple(sorted(vartuple, key=lambda v: v.ptr()))
         self.ptrtuple = tuple(v.ptr() for v in self.vartuple)
-        self.hashval = sum(self.ptrtuple)
+        self.hashval = <Py_ssize_t>hash(self.ptrtuple)
 
     def __getitem__(self, idx):
         return self.vartuple[idx]
 
-    def __hash__(self):
+    def __hash__(self) -> Py_ssize_t:
         return self.hashval
 
-    def __eq__(self, other):
+    def __eq__(self, other: Term):
         return self.ptrtuple == other.ptrtuple
 
     def __len__(self):
@@ -116,6 +129,20 @@ class Term:
     def __repr__(self):
         return 'Term(%s)' % ', '.join([str(v) for v in self.vartuple])
 
+    cpdef double _evaluate(self, Solution sol) except *:
+        cdef double res = 1.0
+        cdef SCIP* scip_ptr = sol.scip
+        cdef SCIP_SOL* sol_ptr = sol.sol
+        cdef int i = 0, n = len(self)
+        cdef Variable var
+
+        for i in range(n):
+            var = <Variable>self.vartuple[i]
+            res *= SCIPgetSolVal(scip_ptr, sol_ptr, var.scip_var)
+            if res == 0:  # early stop
+                return 0.0
+        return res
+
 
 CONST = Term()
 
@@ -157,7 +184,7 @@ def buildGenExprObj(expr):
 ##@details Polynomial expressions of variables with operator overloading. \n
 #See also the @ref ExprDetails "description" in the expr.pxi. 
 cdef class Expr:
-    
+
     def __init__(self, terms=None):
         '''terms is a dict of variables to coefficients.
 
@@ -318,6 +345,20 @@ cdef class Expr:
         else:
             return max(len(v) for v in self.terms)
 
+    cpdef double _evaluate(self, Solution sol) except *:
+        cdef double res = 0
+        cdef Py_ssize_t pos = <Py_ssize_t>0
+        cdef PyObject* key_ptr
+        cdef PyObject* val_ptr
+        cdef Term term
+        cdef double coef
+
+        while PyDict_Next(self.terms, &pos, &key_ptr, &val_ptr):
+            term = <Term>key_ptr
+            coef = <double>(<object>val_ptr)
+            res += coef * term._evaluate(sol)
+        return res
+
 
 cdef class ExprCons:
     '''Constraints with a polynomial expressions and lower/upper bounds.'''
@@ -427,10 +468,10 @@ Operator = Op()
 #
 #See also the @ref ExprDetails "description" in the expr.pxi. 
 cdef class GenExpr:
+
     cdef public _op
     cdef public children
 
-
     def __init__(self): # do we need it
         ''' '''
 
@@ -625,44 +666,88 @@ cdef class SumExpr(GenExpr):
     def __repr__(self):
         return self._op + "(" + str(self.constant) + "," + ",".join(map(lambda child : child.__repr__(), self.children)) + ")"
 
+    cpdef double _evaluate(self, Solution sol) except *:
+        cdef double res = self.constant
+        cdef int i = 0, n = len(self.children)
+        cdef list children = self.children
+        cdef list coefs = self.coefs
+        for i in range(n):
+            res += <double>coefs[i] * (<GenExpr>children[i])._evaluate(sol)
+        return res
+
+
 # Prod Expressions
 cdef class ProdExpr(GenExpr):
+
     cdef public constant
+
     def __init__(self):
         self.constant = 1.0
         self.children = []
         self._op = Operator.prod
+
     def __repr__(self):
         return self._op + "(" + str(self.constant) + "," + ",".join(map(lambda child : child.__repr__(), self.children)) + ")"
 
+    cpdef double _evaluate(self, Solution sol) except *:
+        cdef double res = self.constant
+        cdef list children = self.children
+        cdef int i = 0, n = len(children)
+        for i in range(n):
+            res *= (<GenExpr>children[i])._evaluate(sol)
+            if res == 0:  # early stop
+                return 0.0
+        return res
+
+
 # Var Expressions
 cdef class VarExpr(GenExpr):
+
     cdef public var
+
     def __init__(self, var):
         self.children = [var]
         self._op = Operator.varidx
+
     def __repr__(self):
         return self.children[0].__repr__()
 
+    cpdef double _evaluate(self, Solution sol) except *:
+        return (<Expr>self.children[0])._evaluate(sol)
+
+
 # Pow Expressions
 cdef class PowExpr(GenExpr):
+
     cdef public expo
+
     def __init__(self):
         self.expo = 1.0
         self.children = []
         self._op = Operator.power
+
     def __repr__(self):
         return self._op + "(" + self.children[0].__repr__() + "," + str(self.expo) + ")"
 
+    cpdef double _evaluate(self, Solution sol) except *:
+        return (<GenExpr>self.children[0])._evaluate(sol) ** self.expo
+
+
 # Exp, Log, Sqrt, Sin, Cos Expressions
 cdef class UnaryExpr(GenExpr):
     def __init__(self, op, expr):
         self.children = []
         self.children.append(expr)
         self._op = op
+
     def __repr__(self):
         return self._op + "(" + self.children[0].__repr__() + ")"
 
+    cpdef double _evaluate(self, Solution sol) except *:
+        cdef double res = (<GenExpr>self.children[0])._evaluate(sol)
+        return math.fabs(res) if self._op == "abs" else getattr(math, self._op)(res)
+
+
 # class for constant expressions
 cdef class Constant(GenExpr):
     cdef public number
@@ -673,6 +758,10 @@ cdef class Constant(GenExpr):
     def __repr__(self):
         return str(self.number)
 
+    cpdef double _evaluate(self, Solution sol) except *:
+        return self.number
+
+
 def exp(expr):
     """returns expression with exp-function"""
     if isinstance(expr, MatrixExpr):   
diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi
index 0cfa3a470..ea4436706 100644
--- a/src/pyscipopt/matrix.pxi
+++ b/src/pyscipopt/matrix.pxi
@@ -4,6 +4,7 @@
 """
 from typing import Literal, Optional, Tuple, Union
 import numpy as np
+from numpy.typing import NDArray
 try:
     # NumPy 2.x location
     from numpy.lib.array_utils import normalize_axis_tuple
@@ -12,6 +13,7 @@ except ImportError:
     from numpy.core.numeric import normalize_axis_tuple
 
 cimport numpy as cnp
+from pyscipopt.scip cimport Expr, Solution
 
 cnp.import_array()
 
@@ -142,6 +144,10 @@ class MatrixExpr(np.ndarray):
         return super().__rsub__(other).view(MatrixExpr)
 
 
+    def _evaluate(self, Solution sol) -> NDArray[np.float64]:
+        return _vec_evaluate(self, sol).view(np.ndarray)
+
+
 class MatrixGenExpr(MatrixExpr):
     pass
 
@@ -166,6 +172,9 @@ cdef inline _ensure_array(arg, bool convert_scalar = True):
     return np.array(arg, dtype=object) if convert_scalar else arg
 
 
+_vec_evaluate = np.frompyfunc(lambda expr, sol: expr._evaluate(sol), 2, 1)
+
+
 def _core_dot(cnp.ndarray a, cnp.ndarray b) -> Union[Expr, np.ndarray]:
     """
     Perform matrix multiplication between a N-Demension constant array and a N-Demension
diff --git a/src/pyscipopt/scip.pxd b/src/pyscipopt/scip.pxd
index b9bffc1d6..4743a50a8 100644
--- a/src/pyscipopt/scip.pxd
+++ b/src/pyscipopt/scip.pxd
@@ -2110,6 +2110,8 @@ cdef extern from "tpi/tpi.h":
 cdef class Expr:
     cdef public terms
 
+    cpdef double _evaluate(self, Solution sol)
+
 cdef class Event:
     cdef SCIP_EVENT* event
     # can be used to store problem data
diff --git a/src/pyscipopt/scip.pxi b/src/pyscipopt/scip.pxi
index c4e11ebcb..2895bd380 100644
--- a/src/pyscipopt/scip.pxi
+++ b/src/pyscipopt/scip.pxi
@@ -21,6 +21,7 @@ from dataclasses import dataclass
 from typing import Union
 
 import numpy as np
+from numpy.typing import NDArray
 
 include "expr.pxi"
 include "lp.pxi"
@@ -1099,29 +1100,8 @@ cdef class Solution:
         return sol
 
     def __getitem__(self, expr: Union[Expr, MatrixExpr]):
-        if isinstance(expr, MatrixExpr):
-            result = np.zeros(expr.shape, dtype=np.float64)
-            for idx in np.ndindex(expr.shape):
-                result[idx] = self.__getitem__(expr[idx])
-            return result
-
-        # fast track for Variable
-        cdef SCIP_Real coeff
-        cdef _VarArray wrapper
-        if isinstance(expr, Variable):
-            wrapper = _VarArray(expr)
-            self._checkStage("SCIPgetSolVal")
-            return SCIPgetSolVal(self.scip, self.sol, wrapper.ptr[0])
-        return sum(self._evaluate(term)*coeff for term, coeff in expr.terms.items() if coeff != 0)
-
-    def _evaluate(self, term):
         self._checkStage("SCIPgetSolVal")
-        result = 1
-        cdef _VarArray wrapper
-        wrapper = _VarArray(term.vartuple)
-        for i in range(len(term.vartuple)):
-            result *= SCIPgetSolVal(self.scip, self.sol, wrapper.ptr[i])
-        return result
+        return expr._evaluate(self)
 
     def __setitem__(self, Variable var, value):
         PY_SCIP_CALL(SCIPsetSolVal(self.scip, self.sol, var.scip_var, value))
@@ -10747,7 +10727,11 @@ cdef class Model:
 
         return self.getSolObjVal(self._bestSol, original)
 
-    def getSolVal(self, Solution sol, Expr expr):
+    def getSolVal(
+        self,
+        Solution sol,
+        expr: Union[Expr, GenExpr],
+    ) -> Union[float, NDArray[np.float64]]:
         """
         Retrieve value of given variable or expression in the given solution or in
         the LP/pseudo solution if sol == None
@@ -10767,24 +10751,22 @@ cdef class Model:
         A variable is also an expression.
 
         """
+        if not isinstance(expr, (Expr, GenExpr)):
+            raise TypeError(
+                "Argument 'expr' has incorrect type (expected 'Expr' or 'GenExpr', "
+                f"got {type(expr)})"
+            )
         # no need to create a NULL solution wrapper in case we have a variable
-        cdef _VarArray wrapper
-        if sol == None and isinstance(expr, Variable):
-            wrapper = _VarArray(expr)
-            return SCIPgetSolVal(self._scip, NULL, wrapper.ptr[0])
-        if sol == None:
-            sol = Solution.create(self._scip, NULL)
-        return sol[expr]
+        return (sol or Solution.create(self._scip, NULL))[expr]
 
-    def getVal(self, expr: Union[Expr, MatrixExpr] ):
+    def getVal(self, expr: Union[Expr, GenExpr, MatrixExpr] ):
         """
         Retrieve the value of the given variable or expression in the best known solution.
         Can only be called after solving is completed.
 
         Parameters
         ----------
-        expr : Expr ot MatrixExpr
-            polynomial expression to query the value of
+        expr : Expr, GenExpr or MatrixExpr
 
         Returns
         -------
diff --git a/src/pyscipopt/scip.pyi b/src/pyscipopt/scip.pyi
index 1e703c006..8e39f8d6b 100644
--- a/src/pyscipopt/scip.pyi
+++ b/src/pyscipopt/scip.pyi
@@ -1,5 +1,5 @@
 from dataclasses import dataclass
-from typing import ClassVar
+from typing import TYPE_CHECKING, ClassVar  # noqa: F401
 
 import numpy
 from _typeshed import Incomplete
@@ -2062,7 +2062,6 @@ class Solution:
     data: Incomplete
     def __init__(self, *args: Incomplete, **kwargs: Incomplete) -> None: ...
     def _checkStage(self, method: Incomplete) -> Incomplete: ...
-    def _evaluate(self, term: Incomplete) -> Incomplete: ...
     def getOrigin(self) -> Incomplete: ...
     def retransform(self) -> Incomplete: ...
     def translate(self, target: Incomplete) -> Incomplete: ...
@@ -2122,6 +2121,7 @@ class SumExpr(GenExpr):
     constant: Incomplete
     def __init__(self, *args: Incomplete, **kwargs: Incomplete) -> None: ...
 
+@disjoint_base
 class Term:
     hashval: Incomplete
     ptrtuple: Incomplete
@@ -2138,6 +2138,7 @@ class Term:
     def __lt__(self, other: object) -> bool: ...
     def __ne__(self, other: object) -> bool: ...
 
+@disjoint_base
 class UnaryExpr(GenExpr):
     def __init__(self, *args: Incomplete, **kwargs: Incomplete) -> None: ...
 
diff --git a/tests/test_expr.py b/tests/test_expr.py
index ce79b7cc5..c9135d2fa 100644
--- a/tests/test_expr.py
+++ b/tests/test_expr.py
@@ -1,7 +1,10 @@
+import math
+
 import pytest
 
 from pyscipopt import Model, sqrt, log, exp, sin, cos
-from pyscipopt.scip import Expr, GenExpr, ExprCons, Term, quicksum
+from pyscipopt.scip import Expr, GenExpr, ExprCons, Term
+
 
 @pytest.fixture(scope="module")
 def model():
@@ -188,3 +191,30 @@ def test_rpow_constant_base(model):
 
     with pytest.raises(ValueError):
         c = (-2)**x
+
+
+def test_getVal_with_GenExpr():
+    m = Model()
+    x = m.addVar(lb=1, ub=1, name="x")
+    y = m.addVar(lb=2, ub=2, name="y")
+    z = m.addVar(lb=0, ub=0, name="z")
+    m.optimize()
+
+    # test "Expr({Term(x, y, z): 1.0})"
+    assert m.getVal(z * x * y) == 0
+    # test "Expr({Term(x): 1.0, Term(y): 1.0, Term(): 1.0})"
+    assert m.getVal(x + y + 1) == 4
+    # test "prod(1.0,sum(0.0,prod(1.0,x)),**(sum(0.0,prod(1.0,x)),-1))"
+    assert m.getVal(x / x) == 1
+    # test "prod(1.0,sum(0.0,prod(1.0,y)),**(sum(0.0,prod(1.0,x)),-1))"
+    assert m.getVal(y / x) == 2
+    # test "**(prod(1.0,**(sum(0.0,prod(1.0,x)),-1)),2)"
+    assert m.getVal((1 / x) ** 2) == 1
+    # test "sin(sum(0.0,prod(1.0,x)))"
+    assert round(m.getVal(sin(x)), 6) == round(math.sin(1), 6)
+
+    with pytest.raises(TypeError):
+        m.getVal(1)
+
+    with pytest.raises(ZeroDivisionError):
+        m.getVal(1 / z)
diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py
index 9eeca7eed..b4610a4b7 100644
--- a/tests/test_matrix_variable.py
+++ b/tests/test_matrix_variable.py
@@ -660,3 +660,12 @@ def test_broadcast():
     m.optimize()
 
     assert (m.getVal(x) == np.zeros((2, 3))).all()
+
+
+def test_evaluate():
+    m = Model()
+    x = m.addMatrixVar((1, 1), lb=1, ub=1)
+    m.optimize()
+
+    assert type(m.getVal(x)) is np.ndarray
+    assert m.getVal(x).sum() == 1