diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl
index 962cfe99..0ceb9089 100644
--- a/src/TaylorSeries.jl
+++ b/src/TaylorSeries.jl
@@ -22,7 +22,13 @@ using SparseArrays: SparseMatrixCSC
 using Markdown
 using Requires
 
-import LinearAlgebra: norm, mul!
+using LinearAlgebra: norm, mul!,
+    lu, lu!, LinearAlgebra.lutype, LinearAlgebra.copy_oftype,
+    LinearAlgebra.issuccess
+    # istriu, triu!, istril, tril!, UpperTriangular, LowerTriangular,
+    # LinearAlgebra.inv!, LinearAlgebra.checksquare
+
+import LinearAlgebra: norm, mul!, lu
 
 import Base: ==, +, -, *, /, ^
 
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 1f8bd90c..81dc4a3a 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -560,3 +560,40 @@ function mul!(y::Vector{Taylor1{T}},
 
     return y
 end
+
+
+# Adapted from (Julia v1.2) stdlib/v1.2/LinearAlgebra/src/dense.jl#721-734,
+# licensed under MIT "Expat".
+# Specialize a method of `inv` for Matrix{Taylor1{T}}. Simply, avoid pivoting,
+# since the polynomial field is not an ordered one.
+# function Base.inv(A::StridedMatrix{Taylor1{T}}) where T
+#     checksquare(A)
+#     S = Taylor1{typeof((one(T)*zero(T) + one(T)*zero(T))/one(T))}
+#     AA = convert(AbstractArray{S}, A)
+#     if istriu(AA)
+#         Ai = triu!(parent(inv(UpperTriangular(AA))))
+#     elseif istril(AA)
+#         Ai = tril!(parent(inv(LowerTriangular(AA))))
+#     else
+#         # Do not use pivoting !!
+#         Ai = inv!(lu(AA, Val(false)))
+#         Ai = convert(typeof(parent(Ai)), Ai)
+#     end
+#     return Ai
+# end
+
+# Adapted from (Julia v1.2) stdlib/v1.2/LinearAlgebra/src/lu.jl#240-253
+# and (Julia v1.4.0-dev) stdlib/LinearAlgebra/v1.4/src/lu.jl#270-274,
+# licensed under MIT "Expat".
+# Specialize a method of `lu` for Matrix{Taylor1{T}}, which avoids pivoting,
+# since the polynomial field is not an ordered one.
+# We can't assume an ordered field so we first try without pivoting
+function lu(A::AbstractMatrix{Taylor1{T}}; check::Bool = true) where T
+    S = Taylor1{lutype(T)}
+    F = lu!(copy_oftype(A, S), Val(false); check = false)
+    if issuccess(F)
+        return F
+    else
+        return lu!(copy_oftype(A, S), Val(true); check = check)
+    end
+end
diff --git a/test/onevariable.jl b/test/onevariable.jl
index f1465324..982be1a1 100644
--- a/test/onevariable.jl
+++ b/test/onevariable.jl
@@ -502,14 +502,34 @@ eeuler = Base.MathConstants.e
     @test Taylor1{Float64}(false) == Taylor1([0.0])
     @test Taylor1{Int}(true) == Taylor1([1])
     @test Taylor1{Int}(false) == Taylor1([0])
+end
 
-    # Test use of `inv` for computation of matrix inverse of `Matrix{Taylor1{Float64}}`
-    a = rand(3, 3)
+@testset "Test `inv` for `Matrix{Taylor1{Float64}}``" begin
+    t = Taylor1(5)
+    a = Diagonal(rand(0:10,3)) + rand(3, 3)
+    ainv = inv(a)
     b = Taylor1.(a, 5)
     binv = inv(b)
-    I_t1_5 = Taylor1.(Matrix{Float64}(I, size(b)), 5) # 5x5 Taylor1{Float64} identity matrix, order 5
-    @test norm(b*binv - I_t1_5, Inf) ≤ 1e-12
-    @test norm(binv*b - I_t1_5, Inf) ≤ 1e-12
+    tol = 1.0e-11
+
+    for its = 1:10
+        a .= Diagonal(rand(2:12,3)) + rand(3, 3)
+        ainv .= inv(a)
+        b .= Taylor1.(a, 5)
+        binv .= inv(b)
+        @test norm(binv - ainv, Inf) ≤ tol
+        @test norm(b*binv - I, Inf) ≤ tol
+        @test norm(binv*b - I, Inf) ≤ tol
+        @test norm(triu(b)*inv(UpperTriangular(b)) - I, Inf) ≤ tol
+        @test norm(inv(LowerTriangular(b))*tril(b) - I, Inf) ≤ tol
+
+        b .= b .+ t
+        binv .= inv(b)
+        @test norm(b*binv - I, Inf) ≤ tol
+        @test norm(binv*b - I, Inf) ≤ tol
+        @test norm(triu(b)*inv(triu(b)) - I, Inf) ≤ tol
+        @test norm(inv(tril(b))*tril(b) - I, Inf) ≤ tol
+    end
 end
 
 @testset "Matrix multiplication for Taylor1" begin