diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml
index 44db983..aabc045 100644
--- a/.github/workflows/TagBot.yml
+++ b/.github/workflows/TagBot.yml
@@ -29,5 +29,4 @@ jobs:
       - uses: JuliaRegistries/TagBot@v1
         with:
           token: ${{ secrets.TAGBOT_TOKEN }}
-          # Edit the following line to reflect the actual name of the GitHub Secret containing your private key
           registry: HolyLab/HolyLabRegistry
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 9b4671f..b9aeb12 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "CachedInterpolations"
 uuid = "b9709bfb-d23d-5560-8276-8c35c4b76823"
-version = "1.0.0"
+version = "1.0.1"
 authors = ["Tim Holy <tim.holy@gmail.com>"]
 
 [deps]
@@ -11,7 +11,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Aqua = "0.8"
 Documenter = "1"
 ExplicitImports = "1"
-Interpolations = "0.14 - 1"
+Interpolations = "0.15 - 1"
 StaticArrays = "1"
 Test = "1"
 julia = "1.10"
diff --git a/src/CachedInterpolations.jl b/src/CachedInterpolations.jl
index 66382df..94c75bc 100644
--- a/src/CachedInterpolations.jl
+++ b/src/CachedInterpolations.jl
@@ -131,7 +131,7 @@ julia> itps_c[1](0.0)    # peak via centered coordinates
 0.75
 ```
 """
-function cachedinterpolators(parent::Array{T, M}, N, origin = ntuple(d -> 0, N)) where {T, M}
+function cachedinterpolators(parent::Array{T, M}, N::Integer, origin = ntuple(d -> 0, N)) where {T, M}
     0 <= N <= M || error("N must be between 0 and $M")
     length(origin) == N || throw(DimensionMismatch("length(origin) = $(length(origin)) is inconsistent with $N interpolating dimensions"))
     sz3 = ntuple(d -> d <= N ? 3 : size(parent, d), Val(M))
@@ -172,13 +172,26 @@ end
     return icoefs[wis...]
 end
 
-# FIXME: this function cheats dangerously, because it does _not_
-# update the cache. This is equivalent to the assumption you've called
-# getindex for the current `(x_1, x_2, ...)` location before calling
-# gradient!. If this is not true, you'll get wrong answers.
+"""
+    Interpolations.gradient(itp::CachedInterpolation, ys...)
+
+Return the gradient of `itp` evaluated at coordinates `ys` as an `SVector`.
+The cache is updated automatically, so this can be called without a preceding
+`itp(ys...)`.
+"""
 @inline function Interpolations.gradient(itp::CachedInterpolation{T, N, M, O, K}, ys::Vararg{Number, N}) where {T, N, M, O, K}
-    coefs, tileindex = itp.coefs, itp.tileindex
-    xs = ys .- round.(Int, ys) .+ 2
+    coefs, parent, center, tileindex = itp.coefs, itp.parent, itp.center, itp.tileindex
+    iys = round.(Int, ys)
+    xs = ys .- iys .+ 2
+    newcenter = iys .+ O
+    sz3 = ntuple(d -> 3, Val(N))
+    if newcenter != center
+        offset = CartesianIndex(newcenter .- 2)
+        for i in CartesianIndices(sz3)
+            coefs[i, tileindex] = parent[i + offset, tileindex]
+        end
+        itp.center = newcenter
+    end
     itpinfo = (ntuple(d -> BSpline(Quadratic(InPlace(OnCell()))), Val(N))..., ntuple(d -> NoInterp(), Val(K))...)
     wis = weightedindexes((value_weights, gradient_weights), itpinfo, axes(coefs), (xs..., Tuple(tileindex)...))
     icoefs = InterpGetindex(coefs)