diff --git a/src/derivatives.jl b/src/derivatives.jl
index 3bd78824..3fb10356 100644
--- a/src/derivatives.jl
+++ b/src/derivatives.jl
@@ -1,23 +1,29 @@
 function derivative(A, t, order = 1)
-    order > 2 && throw(DerivativeNotFoundError())
     ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError())
-    order == 1 && return _derivative(A, t, firstindex(A.t) - 1)[1]
-    return ForwardDiff.derivative(t -> _derivative(A, t, firstindex(A.t) - 1)[1], t)
+    iguess = A.idx_prev[]
+
+    return if order == 1
+        val, idx = _derivative(A, t, iguess)
+        A.idx_prev[] = idx
+        val
+    elseif order == 2
+        ForwardDiff.derivative(t -> begin
+                val, idx = _derivative(A, t, iguess)
+                A.idx_prev[] = idx
+                val
+            end, t)
+    else
+        throw(DerivativeNotFoundError())
+    end
 end
 
 function _derivative(A::LinearInterpolation{<:AbstractVector}, t::Number, iguess)
-    idx = searchsortedfirstcorrelated(A.t, t, iguess)
-    idx > length(A.t) ? idx -= 1 : nothing
-    idx -= 1
-    idx == 0 ? idx += 1 : nothing
+    idx = get_idx(A.t, t, iguess; idx_shift = -1, ub_shift = -2, side = :first)
     (A.u[idx + 1] - A.u[idx]) / (A.t[idx + 1] - A.t[idx]), idx
 end
 
 function _derivative(A::LinearInterpolation{<:AbstractMatrix}, t::Number, iguess)
-    idx = searchsortedfirstcorrelated(A.t, t, iguess)
-    idx > length(A.t) ? idx -= 1 : nothing
-    idx -= 1
-    idx == 0 ? idx += 1 : nothing
+    idx = get_idx(A.t, t, iguess; idx_shift = -1, ub_shift = -2, side = :first)
     (@views @. (A.u[:, idx + 1] - A.u[:, idx]) / (A.t[idx + 1] - A.t[idx])), idx
 end
 
@@ -105,17 +111,18 @@ function _derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number)
     der
 end
 
-_derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number, i) = _derivative(A, t), i
-_derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, i) = _derivative(A, t), i
+function _derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number, idx)
+    _derivative(A, t), idx
+end
+function _derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, idx)
+    _derivative(A, t), idx
+end
 
 function _derivative(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess)
-    i = searchsortedfirstcorrelated(A.t, t, iguess)
-    i > length(A.t) ? i -= 1 : nothing
-    i -= 1
-    i == 0 ? i += 1 : nothing
-    j = min(i, length(A.c))  # for smooth derivative at A.t[end]
-    wj = t - A.t[i]
-    (@evalpoly wj A.b[i] 2A.c[j] 3A.d[j]), i
+    idx = get_idx(A.t, t, iguess; idx_shift = -1, side = :first)
+    j = min(idx, length(A.c))  # for smooth derivative at A.t[end]
+    wj = t - A.t[idx]
+    (@evalpoly wj A.b[idx] 2A.c[j] 3A.d[j]), idx
 end
 
 function _derivative(A::ConstantInterpolation{<:AbstractVector}, t::Number)
@@ -130,32 +137,26 @@ end
 
 # QuadraticSpline Interpolation
 function _derivative(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess)
-    idx = searchsortedfirstcorrelated(A.t, t, iguess)
-    idx > length(A.t) ? idx -= 1 : nothing
-    idx == 1 ? idx += 1 : nothing
+    idx = get_idx(A.t, t, iguess; lb = 2, ub_shift = 0, side = :first)
     σ = 1 // 2 * (A.z[idx] - A.z[idx - 1]) / (A.t[idx] - A.t[idx - 1])
     A.z[idx - 1] + 2σ * (t - A.t[idx - 1]), idx
 end
 
 # CubicSpline Interpolation
 function _derivative(A::CubicSpline{<:AbstractVector}, t::Number, iguess)
-    i = searchsortedfirstcorrelated(A.t, t, iguess)
-    i > length(A.t) ? i -= 1 : nothing
-    i -= 1
-    i == 0 ? i += 1 : nothing
-    dI = -3A.z[i] * (A.t[i + 1] - t)^2 / (6A.h[i + 1]) +
-         3A.z[i + 1] * (t - A.t[i])^2 / (6A.h[i + 1])
-    dC = A.u[i + 1] / A.h[i + 1] - A.z[i + 1] * A.h[i + 1] / 6
-    dD = -(A.u[i] / A.h[i + 1] - A.z[i] * A.h[i + 1] / 6)
-    dI + dC + dD, i
+    idx = get_idx(A.t, t, iguess)
+    dI = -3A.z[idx] * (A.t[idx + 1] - t)^2 / (6A.h[idx + 1]) +
+         3A.z[idx + 1] * (t - A.t[idx])^2 / (6A.h[idx + 1])
+    dC = A.u[idx + 1] / A.h[idx + 1] - A.z[idx + 1] * A.h[idx + 1] / 6
+    dD = -(A.u[idx] / A.h[idx + 1] - A.z[idx] * A.h[idx + 1] / 6)
+    dI + dC + dD, idx
 end
 
 function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number, iguess)
     # change t into param [0 1]
     t < A.t[1] && return zero(A.u[1]), 1
     t > A.t[end] && return zero(A.u[end]), lastindex(t)
-    idx = searchsortedlastcorrelated(A.t, t, iguess)
-    idx == length(A.t) ? idx -= 1 : nothing
+    idx = get_idx(A.t, t, iguess)
     n = length(A.t)
     scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
     t_ = A.p[idx] + (t - A.t[idx]) * scale
@@ -176,8 +177,7 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig
     # change t into param [0 1]
     t < A.t[1] && return zero(A.u[1]), 1
     t > A.t[end] && return zero(A.u[end]), lastindex(t)
-    idx = searchsortedlastcorrelated(A.t, t, iguess)
-    idx == length(A.t) ? idx -= 1 : nothing
+    idx = get_idx(A.t, t, iguess)
     scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
     t_ = A.p[idx] + (t - A.t[idx]) * scale
     N = spline_coefficients(A.h, A.d - 1, A.k, t_)
diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index 62352001..86c6be61 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -17,8 +17,9 @@ struct LinearInterpolation{uType, tType, T} <: AbstractInterpolation{T}
     u::uType
     t::tType
     extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
     function LinearInterpolation(u, t, extrapolate)
-        new{typeof(u), typeof(t), eltype(u)}(u, t, extrapolate)
+        new{typeof(u), typeof(t), eltype(u)}(u, t, extrapolate, Ref(1))
     end
 end
 
@@ -48,10 +49,11 @@ struct QuadraticInterpolation{uType, tType, T} <: AbstractInterpolation{T}
     t::tType
     mode::Symbol
     extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
     function QuadraticInterpolation(u, t, mode, extrapolate)
         mode ∈ (:Forward, :Backward) ||
             error("mode should be :Forward or :Backward for QuadraticInterpolation")
-        new{typeof(u), typeof(t), eltype(u)}(u, t, mode, extrapolate)
+        new{typeof(u), typeof(t), eltype(u)}(u, t, mode, extrapolate, Ref(1))
     end
 end
 
@@ -86,6 +88,7 @@ struct LagrangeInterpolation{uType, tType, T, bcacheType} <:
     n::Int
     bcache::bcacheType
     extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
     function LagrangeInterpolation(u, t, n, extrapolate)
         bcache = zeros(eltype(u[1]), n + 1)
         fill!(bcache, NaN)
@@ -93,7 +96,9 @@ struct LagrangeInterpolation{uType, tType, T, bcacheType} <:
             t,
             n,
             bcache,
-            extrapolate)
+            extrapolate,
+            Ref(1)
+        )
     end
 end
 
@@ -128,6 +133,7 @@ struct AkimaInterpolation{uType, tType, bType, cType, dType, T} <:
     c::cType
     d::dType
     extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
     function AkimaInterpolation(u, t, b, c, d, extrapolate)
         new{typeof(u), typeof(t), typeof(b), typeof(c),
             typeof(d), eltype(u)}(u,
@@ -135,7 +141,9 @@ struct AkimaInterpolation{uType, tType, bType, cType, dType, T} <:
             b,
             c,
             d,
-            extrapolate)
+            extrapolate,
+            Ref(1)
+        )
     end
 end
 
@@ -186,8 +194,9 @@ struct ConstantInterpolation{uType, tType, dirType, T} <: AbstractInterpolation{
     t::tType
     dir::Symbol # indicates if value to the $dir should be used for the interpolation
     extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
     function ConstantInterpolation(u, t, dir, extrapolate)
-        new{typeof(u), typeof(t), typeof(dir), eltype(u)}(u, t, dir, extrapolate)
+        new{typeof(u), typeof(t), typeof(dir), eltype(u)}(u, t, dir, extrapolate, Ref(1))
     end
 end
 
@@ -219,6 +228,7 @@ struct QuadraticSpline{uType, tType, tAType, dType, zType, T} <:
     d::dType
     z::zType
     extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
     function QuadraticSpline(u, t, tA, d, z, extrapolate)
         new{typeof(u), typeof(t), typeof(tA),
             typeof(d), typeof(z), eltype(u)}(u,
@@ -226,7 +236,9 @@ struct QuadraticSpline{uType, tType, tAType, dType, zType, T} <:
             tA,
             d,
             z,
-            extrapolate)
+            extrapolate,
+            Ref(1)
+        )
     end
 end
 
@@ -286,12 +298,15 @@ struct CubicSpline{uType, tType, hType, zType, T} <: AbstractInterpolation{T}
     h::hType
     z::zType
     extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
     function CubicSpline(u, t, h, z, extrapolate)
         new{typeof(u), typeof(t), typeof(h), typeof(z), eltype(u)}(u,
             t,
             h,
             z,
-            extrapolate)
+            extrapolate,
+            Ref(1)
+        )
     end
 end
 
@@ -365,6 +380,7 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, T} <:
     pVecType::Symbol
     knotVecType::Symbol
     extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
     function BSplineInterpolation(u,
             t,
             d,
@@ -382,7 +398,9 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, T} <:
             c,
             pVecType,
             knotVecType,
-            extrapolate)
+            extrapolate,
+            Ref(1)
+        )
     end
 end
 
@@ -482,6 +500,7 @@ struct BSplineApprox{uType, tType, pType, kType, cType, T} <:
     pVecType::Symbol
     knotVecType::Symbol
     extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
     function BSplineApprox(u,
             t,
             d,
@@ -501,7 +520,9 @@ struct BSplineApprox{uType, tType, pType, kType, cType, T} <:
             c,
             pVecType,
             knotVecType,
-            extrapolate)
+            extrapolate,
+            Ref(1)
+        )
     end
 end
 
diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl
index 44bd1561..85344c93 100644
--- a/src/interpolation_methods.jl
+++ b/src/interpolation_methods.jl
@@ -1,7 +1,9 @@
-function _interpolate(interp, t)
-    (!interp.extrapolate && (t < interp.t[1] || t > interp.t[end])) &&
+function _interpolate(A, t)
+    ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) &&
         throw(ExtrapolationError())
-    _interpolate(interp, t, firstindex(interp.t) - 1)[1]
+    val, idx_prev = _interpolate(A, t, A.idx_prev[])
+    A.idx_prev[] = idx_prev
+    return val
 end
 
 # Linear Interpolation
@@ -12,7 +14,7 @@ function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, igues
         t1 = t2 = one(eltype(A.t))
         u1 = u2 = one(eltype(A.u))
     else
-        idx = max(1, min(searchsortedlastcorrelated(A.t, t, iguess), length(A.t) - 1))
+        idx = get_idx(A.t, t, iguess)
         t1, t2 = A.t[idx], A.t[idx + 1]
         u1, u2 = A.u[idx], A.u[idx + 1]
     end
@@ -25,7 +27,7 @@ function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, igues
 end
 
 function _interpolate(A::LinearInterpolation{<:AbstractMatrix}, t::Number, iguess)
-    idx = max(1, min(searchsortedlastcorrelated(A.t, t, iguess), length(A.t) - 1))
+    idx = get_idx(A.t, t, iguess)
     θ = (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx])
     return (1 - θ) * A.u[:, idx] + θ * A.u[:, idx + 1], idx
 end
@@ -33,9 +35,7 @@ end
 # Quadratic Interpolation
 _quad_interp_indices(A, t) = _quad_interp_indices(A, t, firstindex(A.t) - 1)
 function _quad_interp_indices(A::QuadraticInterpolation, t::Number, iguess)
-    inner_idx = searchsortedlastcorrelated(A.t, t, iguess)
-    A.mode == :Backward && (inner_idx -= 1)
-    idx = max(1, min(inner_idx, length(A.t) - 2))
+    idx = get_idx(A.t, t, iguess; idx_shift = A.mode == :Backward ? -1 : 0, ub_shift = -2)
     idx, idx + 1, idx + 2
 end
 
@@ -114,61 +114,59 @@ function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number)
     N / D
 end
 
-function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, i)
-    _interpolate(A, t), i
+function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, idx)
+    _interpolate(A, t), idx
 end
 
-function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, i)
-    _interpolate(A, t), i
+function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, idx)
+    _interpolate(A, t), idx
 end
 
 function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess)
-    i = max(1, min(searchsortedlastcorrelated(A.t, t, iguess), length(A.t) - 1))
-    wj = t - A.t[i]
-    (@evalpoly wj A.u[i] A.b[i] A.c[i] A.d[i]), i
+    idx = get_idx(A.t, t, iguess)
+    wj = t - A.t[idx]
+    (@evalpoly wj A.u[idx] A.b[idx] A.c[idx] A.d[idx]), idx
 end
 
 # ConstantInterpolation Interpolation
 function _interpolate(A::ConstantInterpolation{<:AbstractVector}, t::Number, iguess)
     if A.dir === :left
         # :left means that value to the left is used for interpolation
-        i = max(1, searchsortedlastcorrelated(A.t, t, iguess))
-        return A.u[i], i
+        idx = get_idx(A.t, t, iguess; lb = 1, ub_shift = 0)
     else
         # :right means that value to the right is used for interpolation
-        i = min(length(A.t), searchsortedfirstcorrelated(A.t, t, iguess))
-        return A.u[i], i
+        idx = get_idx(A.t, t, iguess; side = :first, lb = 1, ub_shift = 0)
     end
+    A.u[idx], idx
 end
 
 function _interpolate(A::ConstantInterpolation{<:AbstractMatrix}, t::Number, iguess)
     if A.dir === :left
         # :left means that value to the left is used for interpolation
-        i = max(1, searchsortedlastcorrelated(A.t, t, iguess))
-        return A.u[:, i], i
+        idx = get_idx(A.t, t, iguess; lb = 1, ub_shift = 0)
     else
         # :right means that value to the right is used for interpolation
-        i = min(length(A.t), searchsortedfirstcorrelated(A.t, t, iguess))
-        return A.u[:, i], i
+        idx = get_idx(A.t, t, iguess; side = :first, lb = 1, ub_shift = 0)
     end
+    A.u[:, idx], idx
 end
 
 # QuadraticSpline Interpolation
 function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess)
-    i = min(max(2, searchsortedfirstcorrelated(A.t, t, iguess)), length(A.t))
-    Cᵢ = A.u[i - 1]
-    σ = 1 // 2 * (A.z[i] - A.z[i - 1]) / (A.t[i] - A.t[i - 1])
-    return A.z[i - 1] * (t - A.t[i - 1]) + σ * (t - A.t[i - 1])^2 + Cᵢ, i
+    idx = get_idx(A.t, t, iguess; lb = 2, ub_shift = 0, side = :first)
+    Cᵢ = A.u[idx - 1]
+    σ = 1 // 2 * (A.z[idx] - A.z[idx - 1]) / (A.t[idx] - A.t[idx - 1])
+    return A.z[idx - 1] * (t - A.t[idx - 1]) + σ * (t - A.t[idx - 1])^2 + Cᵢ, idx
 end
 
 # CubicSpline Interpolation
 function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess)
-    i = max(1, min(searchsortedlastcorrelated(A.t, t, iguess), length(A.t) - 1))
-    I = A.z[i] * (A.t[i + 1] - t)^3 / (6A.h[i + 1]) +
-        A.z[i + 1] * (t - A.t[i])^3 / (6A.h[i + 1])
-    C = (A.u[i + 1] / A.h[i + 1] - A.z[i + 1] * A.h[i + 1] / 6) * (t - A.t[i])
-    D = (A.u[i] / A.h[i + 1] - A.z[i] * A.h[i + 1] / 6) * (A.t[i + 1] - t)
-    I + C + D, i
+    idx = get_idx(A.t, t, iguess)
+    I = A.z[idx] * (A.t[idx + 1] - t)^3 / (6A.h[idx + 1]) +
+        A.z[idx + 1] * (t - A.t[idx])^3 / (6A.h[idx + 1])
+    C = (A.u[idx + 1] / A.h[idx + 1] - A.z[idx + 1] * A.h[idx + 1] / 6) * (t - A.t[idx])
+    D = (A.u[idx] / A.h[idx + 1] - A.z[idx] * A.h[idx + 1] / 6) * (A.t[idx + 1] - t)
+    I + C + D, idx
 end
 
 # BSpline Curve Interpolation
@@ -178,8 +176,7 @@ function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}},
     t < A.t[1] && return A.u[1], 1
     t > A.t[end] && return A.u[end], lastindex(t)
     # change t into param [0 1]
-    idx = searchsortedlastcorrelated(A.t, t, iguess)
-    idx == length(A.t) ? idx -= 1 : nothing
+    idx = get_idx(A.t, t, iguess)
     t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx])
     n = length(A.t)
     N = spline_coefficients(n, A.d, A.k, t)
@@ -195,10 +192,8 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i
     t < A.t[1] && return A.u[1], 1
     t > A.t[end] && return A.u[end], lastindex(t)
     # change t into param [0 1]
-    idx = searchsortedlastcorrelated(A.t, t, iguess)
-    idx == length(A.t) ? idx -= 1 : nothing
+    idx = get_idx(A.t, t, iguess)
     t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx])
-    n = length(A.t)
     N = spline_coefficients(A.h, A.d, A.k, t)
     ucum = zero(eltype(A.u))
     for i in 1:(A.h)
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index 1b416ab5..bb62e78a 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -41,11 +41,12 @@ function munge_data(u::AbstractVector, t::AbstractVector)
     Tu = Base.nonmissingtype(eltype(u))
     Tt = Base.nonmissingtype(eltype(t))
     @assert length(t) == length(u)
-    non_missing_indices = collect(i
-    for i in 1:length(t)
-    if !ismissing(u[i]) && !ismissing(t[i]))
-    newu = Tu.([u[i] for i in non_missing_indices])
-    newt = Tt.([t[i] for i in non_missing_indices])
+    non_missing_indices = collect(
+        idx for idx in 1:length(t)
+    if !ismissing(u[idx]) && !ismissing(t[idx])
+    )
+    newu = Tu.([u[idx] for idx in non_missing_indices])
+    newt = Tt.([t[idx] for idx in non_missing_indices])
 
     return newu, newt
 end
@@ -54,11 +55,23 @@ function munge_data(U::StridedMatrix, t::AbstractVector)
     TU = Base.nonmissingtype(eltype(U))
     Tt = Base.nonmissingtype(eltype(t))
     @assert length(t) == size(U, 2)
-    non_missing_indices = collect(i
-    for i in 1:length(t)
-    if !any(ismissing, U[:, i]) && !ismissing(t[i]))
-    newUs = [TU.(U[:, i]) for i in non_missing_indices]
-    newt = Tt.([t[i] for i in non_missing_indices])
+    non_missing_indices = collect(
+        idx for idx in 1:length(t)
+    if !any(ismissing, U[:, idx]) && !ismissing(t[idx])
+    )
+    newUs = [TU.(U[:, idx]) for idx in non_missing_indices]
+    newt = Tt.([t[idx] for idx in non_missing_indices])
 
     return hcat(newUs...), newt
 end
+
+function get_idx(tvec, t, iguess; lb = 1, ub_shift = -1, idx_shift = 0, side = :last)
+    ub = length(tvec) + ub_shift
+    return if side == :last
+        clamp(searchsortedlastcorrelated(tvec, t, iguess) + idx_shift, lb, ub)
+    elseif side == :first
+        clamp(searchsortedfirstcorrelated(tvec, t, iguess) + idx_shift, lb, ub)
+    else
+        error("side must be :first or :last")
+    end
+end
diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl
index 484c0330..d3abbdd2 100644
--- a/test/derivative_tests.jl
+++ b/test/derivative_tests.jl
@@ -1,4 +1,5 @@
 using DataInterpolations, Test
+using FindFirstFunctions: searchsortedfirstcorrelated
 using FiniteDifferences
 using DataInterpolations: derivative
 using Symbolics
@@ -35,6 +36,11 @@ function test_derivatives(method, u, t; args = [], kwargs = [], name::String)
             adiff2 = derivative(func, _t, 2)
             @test isapprox(fdiff, adiff, atol = 1e-8)
             @test isapprox(fdiff2, adiff2, atol = 1e-8)
+            # Cached index
+            if hasproperty(func, :idx_prev)
+                @test abs(func.idx_prev[] -
+                          searchsortedfirstcorrelated(func.t, _t, func.idx_prev[])) <= 1
+            end
         end
 
         # t = t0
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index a11d6ea4..a7c4eb2e 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -1,7 +1,17 @@
 using DataInterpolations
+using FindFirstFunctions: searchsortedfirstcorrelated
 using StableRNGs
 using Optim, ForwardDiff
 
+function test_cached_index(A)
+    for t in range(first(A.t), last(A.t); length = 2 * length(A.t) - 1)
+        A(t)
+        idx = searchsortedfirstcorrelated(A.t, t, A.idx_prev[])
+        @test abs(A.idx_prev[] -
+                  searchsortedfirstcorrelated(A.t, t, A.idx_prev[])) <= 2
+    end
+end
+
 @testset "Linear Interpolation" begin
     for t in (1.0:10.0, 1.0collect(1:10))
         u = 2.0collect(1:10)
@@ -45,6 +55,7 @@ using Optim, ForwardDiff
     @test A(0) == [-2.0 0.0; -3.0 0.0; -4.0 0.0]
     @test A(3) == [4.0 6.0; 6.0 9.0; 8.0 12.0]
     @test A(5) == [8.0 10.0; 12.0 15.0; 16.0 20.0]
+    test_cached_index(A)
 
     # with NaNs (#113)
     u = [NaN, 1.0, 2.0, 3.0]
@@ -168,6 +179,7 @@ end
     @test A(2.5) == 6.25
     @test A(3.5) == 12.25
     @test A(5.0) == 25
+    test_cached_index(A)
 
     # backward-looking interpolation
     u = [1.0, 4.0, 9.0, 16.0]
@@ -182,6 +194,7 @@ end
     @test A(2.5) == 6.25
     @test A(3.5) == 12.25
     @test A(5.0) == 25
+    test_cached_index(A)
 
     # Test both forward and backward-looking quadratic interpolation
     u = [1.0, 4.5, 6.0, 2.0]
@@ -204,6 +217,9 @@ end
     @test A_f(3.5) ≈ l₂ * u[2] + l₁ * u[3] + l₀ * u[4]
     @test A_b(3.5) ≈ l₂ * u[2] + l₁ * u[3] + l₀ * u[4]
 
+    test_cached_index(A_f)
+    test_cached_index(A_b)
+
     # Matrix interpolation test
     u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0]
     A = QuadraticInterpolation(u, t; extrapolate = true)
@@ -321,6 +337,7 @@ end
     @test A(8.6) ≈ 4.1796554159017080820603951
     @test A(9.9) ≈ 3.4110386597938129327189927
     @test A(10.0) ≈ 3.0
+    test_cached_index(A)
 
     # Test extrapolation
     A = AkimaInterpolation(u, t; extrapolate = true)
@@ -345,6 +362,7 @@ end
         @test A(3.5) == u[1]
         @test A(4.0) == u[1]
         @test A(4.5) == u[1]
+        test_cached_index(A)
 
         A = ConstantInterpolation(u, t; extrapolate = true) # dir=:left is default
         @test A(0.5) == u[1]
@@ -356,6 +374,7 @@ end
         @test A(3.5) == u[3]
         @test A(4.0) == u[1]
         @test A(4.5) == u[1]
+        test_cached_index(A)
     end
 
     @testset "Matrix case" for u in [
@@ -372,6 +391,7 @@ end
         @test A(3.5) == u[:, 1]
         @test A(4.0) == u[:, 1]
         @test A(4.5) == u[:, 1]
+        test_cached_index(A)
 
         A = ConstantInterpolation(u, t; extrapolate = true) # dir=:left is default
         @test A(0.5) == u[:, 1]
@@ -383,6 +403,7 @@ end
         @test A(3.5) == u[:, 3]
         @test A(4.0) == u[:, 1]
         @test A(4.5) == u[:, 1]
+        test_cached_index(A)
     end
 
     @testset "Vector of Vectors case" for u in [
@@ -398,6 +419,7 @@ end
         @test A(3.5) == u[4]
         @test A(4.0) == u[4]
         @test A(4.5) == u[4]
+        test_cached_index(A)
 
         A = ConstantInterpolation(u, t; extrapolate = true) # dir=:left is default
         @test A(0.5) == u[1]
@@ -409,6 +431,7 @@ end
         @test A(3.5) == u[3]
         @test A(4.0) == u[4]
         @test A(4.5) == u[4]
+        test_cached_index(A)
     end
 
     @testset "Vector of Matrices case" for u in [
@@ -424,6 +447,7 @@ end
         @test A(3.5) == u[4]
         @test A(4.0) == u[4]
         @test A(4.5) == u[4]
+        test_cached_index(A)
 
         A = ConstantInterpolation(u, t; extrapolate = true) # dir=:left is default
         @test A(0.5) == u[1]
@@ -435,6 +459,7 @@ end
         @test A(3.5) == u[3]
         @test A(4.0) == u[4]
         @test A(4.5) == u[4]
+        test_cached_index(A)
     end
 
     # Test extrapolation
@@ -464,6 +489,7 @@ end
     @test A(-0.5) == P₁(-0.5)
     @test A(0.7) == P₂(0.7)
     @test A(2.0) == P₂(2.0)
+    test_cached_index(A)
 
     u_ = [0.0, 1.0, 3.0]' .* ones(4)
     u = [u_[:, i] for i in 1:size(u_, 2)]
@@ -496,6 +522,7 @@ end
     t = [-1.0, 0.0, 1.0]
 
     A = CubicSpline(u, t; extrapolate = true)
+    test_cached_index(A)
 
     # Solution
     P₁ = x -> 1 + 1.5x + 0.75 * x^2 + 0.25 * x^3 # for x ∈ [-1.0, 0.0]
@@ -552,6 +579,7 @@ end
         @test [A(25.0), A(80.0)] == [13.454197730061425, 10.305633616059845]
         @test [A(190.0), A(225.0)] == [14.07428439395079, 11.057784141519251]
         @test [A(t[1]), A(t[end])] == [u[1], u[end]]
+        test_cached_index(A)
 
         # Test extrapolation
         A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform; extrapolate = true)
@@ -582,6 +610,7 @@ end
         @test [A(25.0), A(80.0)] ≈ [12.979802931218234, 10.914310609953178]
         @test [A(190.0), A(225.0)] ≈ [13.851245975109263, 12.963685868886575]
         @test [A(t[1]), A(t[end])] ≈ [u[1], u[end]]
+        test_cached_index(A)
 
         # Test extrapolation
         A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform; extrapolate = true)