diff --git a/docs/src/free_associative_algebra.md b/docs/src/free_associative_algebra.md
index cca5add4df..840a850d9a 100644
--- a/docs/src/free_associative_algebra.md
+++ b/docs/src/free_associative_algebra.md
@@ -167,6 +167,10 @@ julia> exponent_word(-x^2*y^2*z^2*x + z*y, 1)
  1
 ```
 
+```@docs
+evaluate(a::AbstractAlgebra.FreeAssAlgElem{T}, vals::Vector{U}) where {T <: RingElement, U <: NCRingElem}
+```
+
 ### Iterators
 
 The following iterators are provided for elements of a free associative algebra,
diff --git a/src/FreeAssAlgebra.jl b/src/FreeAssAlgebra.jl
index 9f3ab1b5e6..fc9306b971 100644
--- a/src/FreeAssAlgebra.jl
+++ b/src/FreeAssAlgebra.jl
@@ -149,6 +149,66 @@ function Base.hash(x::FreeAssAlgElem{T}, h::UInt) where T <: RingElement
    return b
 end
 
+###############################################################################
+#
+#   Evaluation
+#
+###############################################################################
+
+@doc raw"""
+    evaluate(a::AbstractAlgebra.FreeAssAlgElem{T}, vals::Vector{U}) where {T <: RingElement, U <: NCRingElem}
+
+Evaluate `a` by substituting in the array of values for each of the variables.
+The evaluation will succeed if multiplication is defined between elements of
+the coefficient ring of `a` and elements of `vals`.
+
+The syntax `a(vals...)` is also supported.
+
+# Examples
+
+```jldoctest; setup = :(using AbstractAlgebra)
+julia> R, (x, y) = free_associative_algebra(ZZ, ["x", "y"]);
+
+julia> f = x*y - y*x
+x*y - y*x
+
+julia> S = MatrixAlgebra(ZZ, 2);
+
+julia> m1 = S([1 2; 3 4])
+[1   2]
+[3   4]
+
+julia> m2 = S([0 1; 1 0])
+[0   1]
+[1   0]
+
+julia> evaluate(f, [m1, m2])
+[-1   -3]
+[ 3    1]
+
+julia> m1*m2 - m2*m1 == evaluate(f, [m1, m2])
+true
+
+julia> m1*m2 - m2*m1 == f(m1, m2)
+true
+```
+"""
+function evaluate(a::AbstractAlgebra.FreeAssAlgElem{T}, vals::Vector{U}) where {T <: RingElement, U <: NCRingElem}
+   length(vals) != nvars(parent(a)) && error("Number of variables does not match number of values")
+   R = base_ring(parent(a))
+   S = parent(one(R)*one(parent(vals[1])))
+   r = zero(S)
+   o = one(S)
+   for (c, v) in zip(coefficients(a), exponent_words(a))
+      r = addeq!(r, c*prod((vals[i] for i in v), init = o))
+   end
+   return r
+end
+
+function (a::AbstractAlgebra.FreeAssAlgElem{T})(vals::U...) where {T <: RingElement, U <: NCRingElem}
+   return evaluate(a, [vals...])
+end
+
 ###############################################################################
 #
 #   Random elements