diff --git a/Project.toml b/Project.toml
index 52e8c6dd..94b67700 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "FiniteDifferences"
 uuid = "26cc04aa-876d-5657-8c51-4c34ba976000"
-version = "0.9.5"
+version = "0.9.6"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/methods.jl b/src/methods.jl
index 41485561..dc6838c9 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -156,15 +156,20 @@ function _check_p_q(p::Integer, q::Integer)
 end
 
 # Compute coefficients for the method
+
+const _COEFFS_CACHE = Dict{Tuple{AbstractVector{<:Real}, Integer, Integer}, Vector{Float64}}()
 function _coefs(grid::AbstractVector{<:Real}, p::Integer, q::Integer)
-    # For high precision on the \ we use Rational, and to prevent overflows we use Int128
-    # At the end we go to Float64 for fast floating point math (rather than rational math)
-    C = [Rational{Int128}(g^i) for i in 0:(p - 1), g in grid]
-    x = zeros(Rational{Int128}, p)
-    x[q + 1] = factorial(q)
-    return Float64.(C \ x)
+    return get!(_COEFFS_CACHE, (grid, p, q)) do
+        # For high precision on the \ we use Rational, and to prevent overflows we use Int128
+        # At the end we go to Float64 for fast floating point math (rather than rational math)
+        C = [Rational{Int128}(g^i) for i in 0:(p - 1), g in grid]
+        x = zeros(Rational{Int128}, p)
+        x[q + 1] = factorial(q)
+        return Float64.(C \ x)
+    end
 end
 
+
 # Estimate the bound on the function value and its derivatives at a point
 _estimate_bound(x, cond) = cond * maximum(abs, x) + TINY