diff --git a/Project.toml b/Project.toml
index a7f2851..278b7c8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8"
 
 [compat]
 julia = "1"
diff --git a/src/IntervalMatrices.jl b/src/IntervalMatrices.jl
index 36da1eb..04ead67 100644
--- a/src/IntervalMatrices.jl
+++ b/src/IntervalMatrices.jl
@@ -6,6 +6,8 @@ using LinearAlgebra
 import Random
 using Random: AbstractRNG, GLOBAL_RNG, seed!
 
+import SymEngine
+
 using Reexport
 @reexport using IntervalArithmetic
 
@@ -18,6 +20,7 @@ include("matrix.jl")
 # Arithmetic for interval matrices
 # =================================
 include("arithmetic.jl")
+include("power.jl")
 
 # =======================================================
 # Methods to handle the exponential of an interval matrix
diff --git a/src/power.jl b/src/power.jl
new file mode 100644
index 0000000..b87953e
--- /dev/null
+++ b/src/power.jl
@@ -0,0 +1,31 @@
+import Base: ^
+
+function ^(M::IntervalMatrix{T}, k::Integer) where {T}
+    n = checksquare(M)
+
+    # create symbolic matrix
+    Ms = Matrix{SymEngine.Basic}(undef, n, n)
+    dict = Dict{SymEngine.Basic, Interval{T}}()
+    dict2 = Dict{SymEngine.Basic, T}()
+    for i in 1:n, j in 1:n
+        var = SymEngine.symbols("M_$(i)_$(j)")
+        Ms[i, j] = var
+        dict[var] = M[i, j]
+        dict2[var] = mid(M[i, j])
+    end
+
+    # compute symbolic matrix power
+    Msᵏ = Ms^k
+
+    # substitute intervals in symbolic result
+    Mᵏ = similar(M)
+    for i in 1:n, j in 1:n
+        expr = Msᵏ[i, j]
+#         for (k, v) in dict
+#             expr = SymEngine.subs(expr, k, mid(v))
+#         end
+        expr = SymEngine.subs(expr, dict2...)
+        Mᵏ[i, j] = SymEngine.lambdify(expr)
+    end
+    return Mᵏ
+end