diff --git a/Project.toml b/Project.toml
index 531144a1c..d5c5bfc21 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,7 +19,6 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
 TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
-Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
@@ -36,7 +35,6 @@ Reexport = "0.2, 1.0"
 StructArrays = "0.4,0.5, 0.6"
 TaylorIntegration = "0.8, 0.9"
 TaylorSeries = "0.10, 0.11, 0.12"
-Tullio = "0.3"
 UnPack = "1"
 julia = "1.6,1.7"
 
diff --git a/src/ProbNumDiffEq.jl b/src/ProbNumDiffEq.jl
index f4385849a..0334878f6 100644
--- a/src/ProbNumDiffEq.jl
+++ b/src/ProbNumDiffEq.jl
@@ -32,7 +32,6 @@ using UnPack
 using RecipesBase
 using RecursiveArrayTools
 using ForwardDiff
-using Tullio
 using Octavian
 # By default use mul!
 _matmul!(C, A, B) = mul!(C, A, B)
@@ -48,8 +47,6 @@ _matmul!(
     A::Diagonal{T},
     B::AbstractMatrix{T},
 ) where {T<:LinearAlgebra.BlasFloat} = (C .= A.diag .* B)
-_matmul!(C::Diagonal, A::AbstractMatrix, B::AbstractMatrix) =
-    @tullio C[i, i] = A[i, j] * B[j, i]
 _matmul!(
     C::AbstractMatrix{T},
     A::AbstractVecOrMat{T},