diff --git a/CHANGELOG.md b/CHANGELOG.md
index f82f815f..2870eec6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,7 @@ All notable changes to this project will be documented in this file. The format
 - Add an optional relative-residuals argument to `ConstitutiveMaterial.optimize(relative=False)`.
 - Add a class-decorator `constitutive_material(Msterial, name=None)`.
 - Add the MORPH material formulation implemented by the concept of representative directions `morph_representative_directions(p)` to be used in `Hyperelastic()`.
+- Add the MORPH material formulation `morph(p)` to be used in `Hyperelastic()`.
 
 ### Changed
 - Recfactor the `constitution` module.
diff --git a/docs/felupe/constitution.rst b/docs/felupe/constitution.rst
index e0121bc3..ee352200 100644
--- a/docs/felupe/constitution.rst
+++ b/docs/felupe/constitution.rst
@@ -78,6 +78,7 @@ There are many different pre-defined constitutive material formulations availabl
    finite_strain_viscoelastic
    miehe_goektepe_lulei
    mooney_rivlin
+   morph
    morph_representative_directions
    neo_hooke
    ogden
@@ -229,6 +230,8 @@ There are many different pre-defined constitutive material formulations availabl
 
 .. autofunction:: felupe.mooney_rivlin
 
+.. autofunction:: felupe.morph
+
 .. autofunction:: felupe.morph_representative_directions
 
 .. autofunction:: felupe.neo_hooke
diff --git a/pyproject.toml b/pyproject.toml
index 911cecfe..5c229879 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -52,7 +52,7 @@ dynamic = ["version"]
 dependencies = [
   "numpy",
   "scipy",
-  "tensortrax",
+  "tensortrax>=0.21.1",
 ]
 
 [project.optional-dependencies]
diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py
index 64de5068..b0017335 100644
--- a/src/felupe/__init__.py
+++ b/src/felupe/__init__.py
@@ -46,6 +46,7 @@
     linear_elastic_plastic_isotropic_hardening,
     miehe_goektepe_lulei,
     mooney_rivlin,
+    morph,
     morph_representative_directions,
     neo_hooke,
     ogden,
@@ -188,6 +189,7 @@
     "isochoric_volumetric_split",
     "miehe_goektepe_lulei",
     "mooney_rivlin",
+    "morph",
     "morph_representative_directions",
     "neo_hooke",
     "ogden",
diff --git a/src/felupe/constitution/__init__.py b/src/felupe/constitution/__init__.py
index d4b9ac43..edcd0f7d 100644
--- a/src/felupe/constitution/__init__.py
+++ b/src/felupe/constitution/__init__.py
@@ -19,6 +19,7 @@
     isochoric_volumetric_split,
     miehe_goektepe_lulei,
     mooney_rivlin,
+    morph,
     morph_representative_directions,
     neo_hooke,
     ogden,
@@ -52,6 +53,7 @@
     "isochoric_volumetric_split",
     "miehe_goektepe_lulei",
     "mooney_rivlin",
+    "morph",
     "morph_representative_directions",
     "neo_hooke",
     "ogden",
diff --git a/src/felupe/constitution/hyperelasticity/models/__init__.py b/src/felupe/constitution/hyperelasticity/models/__init__.py
index b76f5ba9..0265085d 100644
--- a/src/felupe/constitution/hyperelasticity/models/__init__.py
+++ b/src/felupe/constitution/hyperelasticity/models/__init__.py
@@ -16,7 +16,7 @@
 from ._helpers import isochoric_volumetric_split
 from ._miehe_goektepe_lulei import miehe_goektepe_lulei
 from ._mooney_rivlin import mooney_rivlin
-from ._morph import morph_representative_directions
+from ._morph import morph, morph_representative_directions
 from ._neo_hooke import neo_hooke
 from ._ogden import ogden
 from ._ogden_roxburgh import ogden_roxburgh
@@ -33,6 +33,7 @@
     "isochoric_volumetric_split",
     "miehe_goektepe_lulei",
     "mooney_rivlin",
+    "morph",
     "morph_representative_directions",
     "neo_hooke",
     "ogden",
diff --git a/src/felupe/constitution/hyperelasticity/models/_morph.py b/src/felupe/constitution/hyperelasticity/models/_morph.py
index 06ab641c..c41cc7f1 100644
--- a/src/felupe/constitution/hyperelasticity/models/_morph.py
+++ b/src/felupe/constitution/hyperelasticity/models/_morph.py
@@ -16,8 +16,118 @@
 along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
 """
 from .microsphere import affine_stretch_statevars
-from tensortrax.math import array, abs as tensor_abs, maximum, sqrt, exp
-from tensortrax.math.special import try_stack
+from tensortrax.math import (
+    array,
+    abs as tensor_abs,
+    maximum,
+    sqrt,
+    exp,
+    if_else,
+    real_to_dual,
+)
+from tensortrax.math.special import try_stack, from_triu_1d, triu_1d, sym, dev, ddot
+from tensortrax.math.linalg import det, inv, eigvalsh, expm
+
+
+def morph(C, statevars, p):
+    """Strain energy function of the
+    `MORPH <https://doi.org/10.1016/s0749-6419(02)00091-8>`_ model formulation [1]_.
+
+    Parameters
+    ----------
+    C : tensortrax.Tensor
+        Right Cauchy-Green deformation tensor.
+    statevars : array
+        Vector of stacked state variables (CTS, C, SA).
+    p : list of float
+        A list which contains the 8 material parameters.
+
+    Examples
+    --------
+    ..  pyvista-plot::
+        :context:
+
+        >>> import felupe as fem
+        >>>
+        >>> umat = fem.Hyperelastic(
+        ...     fem.morph,
+        ...     p=[0.039, 0.371, 0.174, 2.41, 0.0094, 6.84, 5.65, 0.244],
+        ...     nstatevars=13,
+        ... )
+        >>> ax = umat.plot(
+        ...    incompressible=True,
+        ...    ux=fem.math.linsteps(
+        ...        [1, 2, 1, 2.75, 1, 3.5, 1, 4.2, 1, 4.8, 1, 4.8, 1],
+        ...        num=50,
+        ...    ),
+        ...    ps=None,
+        ...    bx=None,
+        ... )
+
+    ..  pyvista-plot::
+        :include-source: False
+        :context:
+        :force_static:
+
+        >>> import pyvista as pv
+        >>>
+        >>> fig = ax.get_figure()
+        >>> chart = pv.ChartMPL(fig)
+        >>> chart.show()
+
+    References
+    ----------
+    .. [1] D. Besdo and J. Ihlemann, "A phenomenological constitutive model for
+       rubberlike materials and its numerical applications", International Journal
+       of Plasticity, vol. 19, no. 7. Elsevier BV, pp. 1019–1036, Jul. 2003. doi:
+       `10.1016/s0749-6419(02)00091-8 <https://doi.org/10.1016/s0749-6419(02)00091-8>`_.
+    """
+
+    # extract old state variables
+    CTSn = array(statevars[0], like=C[0, 0])
+    Cn = from_triu_1d(statevars[1:7], like=C)
+    SAn = from_triu_1d(statevars[7:], like=C)
+
+    # distortional part of right Cauchy-Green deformation tensor
+    I3 = det(C)
+    CG = C * I3 ** (-1 / 3)
+
+    # inverse of and incremental right Cauchy-Green deformation tensor
+    invC = inv(C)
+    dC = C - Cn
+
+    # eigenvalues of right Cauchy-Green deformation tensor (sorted in ascending order)
+    λCG = eigvalsh(CG)
+
+    # Tresca invariant of distortional part of right Cauchy-Green deformation tensor
+    CTG = λCG[-1] - λCG[0]
+
+    # maximum Tresca invariant in load history
+    CTS = maximum(CTG, CTSn)
+
+    def sigmoid(x):
+        "Algebraic sigmoid function."
+        return 1 / sqrt(1 + x**2)
+
+    # material parameters
+    α = p[0] + p[1] * sigmoid(p[2] * CTS)
+    β = p[3] * sigmoid(p[2] * CTS)
+    γ = p[4] * CTS * (1 - sigmoid(CTS / p[5]))
+
+    LG = sym(dev(invC @ dC)) @ CG
+    λLG = eigvalsh(LG)
+    LTG = λLG[-1] - λLG[0]
+    LG_LTG = if_else(LTG > 0, LG / LTG, LG)
+
+    # limiting stresses "L" and additional stresses "A"
+    SL = (γ * expm(p[6] * LG_LTG * CTG / CTS) + p[7] * LG_LTG) @ invC
+    SA = (SAn + β * LTG * SL) / (1 + β * LTG)
+
+    # second Piola-Kirchhoff stress tensor
+    dψdC = α * dev(CG) @ invC + dev(SA @ C) @ invC / 2
+    statevars_new = try_stack([[CTS], triu_1d(C), triu_1d(SA)], fallback=statevars)
+
+    return real_to_dual(dψdC, C, mul=ddot), statevars_new
 
 
 def morph_representative_directions(C, statevars, p, ε=1e-8):
diff --git a/tests/test_constitution.py b/tests/test_constitution.py
index 790ba03f..4f43f70d 100644
--- a/tests/test_constitution.py
+++ b/tests/test_constitution.py
@@ -316,6 +316,13 @@ def test_umat_hyperelastic_statevars():
             ),
             True,
         ),
+        (
+            fem.constitution.morph,
+            dict(
+                p=[0.039, 0.371, 0.174, 2.41, 0.0094, 6.84, 5.65, 0.244], nstatevars=13
+            ),
+            True,
+        ),
     ]:
         umat = fem.Hyperelastic(model, **kwargs)