diff --git a/CHANGELOG.md b/CHANGELOG.md
index 52df46228..65f209b36 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,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 supported GenExpr type
 - Fixed lotsizing_lazy example
 ### Changed
 - changed default value of enablepricing flag to True
diff --git a/src/pyscipopt/expr.pxi b/src/pyscipopt/expr.pxi
index f0c406fcb..84a3d6d99 100644
--- a/src/pyscipopt/expr.pxi
+++ b/src/pyscipopt/expr.pxi
@@ -42,8 +42,16 @@
 # 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
+
 include "matrix.pxi"
 
+if TYPE_CHECKING:
+    double = float
+
 
 def _is_number(e):
     try:
@@ -87,15 +95,17 @@ 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 int 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 = hash(self.ptrtuple)
 
     def __getitem__(self, idx):
         return self.vartuple[idx]
@@ -103,7 +113,7 @@ class Term:
     def __hash__(self):
         return self.hashval
 
-    def __eq__(self, other):
+    def __eq__(self, other: Term):
         return self.ptrtuple == other.ptrtuple
 
     def __len__(self):
@@ -116,6 +126,13 @@ class Term:
     def __repr__(self):
         return 'Term(%s)' % ', '.join([str(v) for v in self.vartuple])
 
+    cpdef double _evaluate(self, Solution sol):
+        cdef double res = 1
+        cdef Variable i
+        for i in self.vartuple:
+            res *= SCIPgetSolVal(sol.scip, sol.sol, i.scip_var)
+        return res
+
 
 CONST = Term()
 
@@ -157,7 +174,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 +335,14 @@ cdef class Expr:
         else:
             return max(len(v) for v in self.terms)
 
+    cpdef double _evaluate(self, Solution sol):
+        cdef double res = 0
+        cdef double coef
+        cdef Term term
+        for term, coef in self.terms.items():
+            res += coef * term._evaluate(sol)
+        return res
+
 
 cdef class ExprCons:
     '''Constraints with a polynomial expressions and lower/upper bounds.'''
@@ -427,10 +452,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 +650,83 @@ 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):
+        cdef double res = self.constant
+        cdef GenExpr child
+        cdef double coef
+        for child, coef in zip(self.children, self.coefs):
+            res += coef * child._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):
+        cdef double res = self.constant
+        cdef GenExpr child
+        for child in self.children:
+            res *= child._evaluate(sol)
+        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):
+        return 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):
+        return 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):
+        return getattr(math, self._op)(self.children[0]._evaluate(sol))
+
+
 # class for constant expressions
 cdef class Constant(GenExpr):
     cdef public number
@@ -673,6 +737,10 @@ cdef class Constant(GenExpr):
     def __repr__(self):
         return str(self.number)
 
+    cpdef double _evaluate(self, Solution sol):
+        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 1a6a09cf3..febdeaf84 100644
--- a/src/pyscipopt/matrix.pxi
+++ b/src/pyscipopt/matrix.pxi
@@ -5,6 +5,7 @@
 
 from typing import 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,8 @@ except ImportError:
     # Fallback for NumPy 1.x
     from numpy.core.numeric import normalize_axis_tuple
 
+from pyscipopt.scip cimport Expr, Solution
+
 
 def _is_number(e):
     try:
@@ -49,6 +52,9 @@ def _matrixexpr_richcmp(self, other, op):
     return res.view(MatrixExprCons)
 
 
+_evaluate = np.frompyfunc(lambda expr, sol: expr._evaluate(sol), 2, 1)
+
+
 class MatrixExpr(np.ndarray):
 
     def sum(
@@ -148,9 +154,14 @@ class MatrixExpr(np.ndarray):
     def __matmul__(self, other):
         return super().__matmul__(other).view(MatrixExpr)
 
+    def _evaluate(self, Solution sol) -> NDArray[np.float64]:
+        return _evaluate(self, sol).view(np.ndarray)
+
+
 class MatrixGenExpr(MatrixExpr):
     pass
 
+
 class MatrixExprCons(np.ndarray):
 
     def __le__(self, other: Union[float, int, np.ndarray]) -> MatrixExprCons:
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 da23028f9..0ab226915 100644
--- a/src/pyscipopt/scip.pxi
+++ b/src/pyscipopt/scip.pxi
@@ -1099,29 +1099,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 +10726,7 @@ 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]) -> float:
         """
         Retrieve value of given variable or expression in the given solution or in
         the LP/pseudo solution if sol == None
@@ -10767,14 +10746,13 @@ 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] ):
         """
diff --git a/src/pyscipopt/scip.pyi b/src/pyscipopt/scip.pyi
index 831dd02ed..f373e6ee0 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
@@ -2065,7 +2065,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: ...
@@ -2125,6 +2124,7 @@ class SumExpr(GenExpr):
     constant: Incomplete
     def __init__(self, *args: Incomplete, **kwargs: Incomplete) -> None: ...
 
+@disjoint_base
 class Term:
     hashval: Incomplete
     ptrtuple: Incomplete
@@ -2141,6 +2141,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..e6b11cdaf 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,21 @@ def test_rpow_constant_base(model):
 
     with pytest.raises(ValueError):
         c = (-2)**x
+
+
+def test_evaluate():
+    m = Model()
+    x = m.addVar(lb=1, ub=1)
+    m.optimize()
+
+    # test "Expr({Term(x1): 1.0, Term(): 1.0})"
+    assert m.getVal(x + 1) == 2
+    # test "prod(1.0,sum(0.0,prod(1.0,x1)),**(sum(0.0,prod(1.0,x1)),-1))"
+    assert m.getVal(x / x) == 1
+    # test "**(prod(1.0,**(sum(0.0,prod(1.0,x1)),-1)),2)"
+    assert m.getVal((1 / x) ** 2) == 1
+    # test "sin(sum(0.0,prod(1.0,x1)))"
+    assert round(m.getVal(sin(x)), 6) == round(math.sin(1), 6)
+
+    with pytest.raises(TypeError):
+        m.getVal(1)
diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py
index e4758f077..c232abcb6 100644
--- a/tests/test_matrix_variable.py
+++ b/tests/test_matrix_variable.py
@@ -604,3 +604,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