diff --git a/Project.toml b/Project.toml
index 544f1654..b06d2adb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+ReadOnlyArrays = "988b38a3-91fc-5605-94a2-ee2116b3bd83"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
@@ -32,6 +33,7 @@ LinearAlgebra = "1.10"
 Optim = "1.6"
 PrettyTables = "2"
 QuadGK = "2.9.1"
+ReadOnlyArrays = "0.2.0"
 RecipesBase = "1.3"
 Reexport = "1"
 RegularizationTools = "0.6"
diff --git a/docs/src/interface.md b/docs/src/interface.md
index 3e7df7e1..ca5e9819 100644
--- a/docs/src/interface.md
+++ b/docs/src/interface.md
@@ -17,7 +17,7 @@ t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3]
 
 All interpolation methods return an object from which we can compute the value of the dependent variable at any time point.
 
-We will use the `CubicSpline` method for demonstration, but the API is the same for all the methods. We can also pass the `extrapolate=true` keyword if we want to allow the interpolation to go beyond the range of the timepoints. The default value is `extrapolate=false`.
+We will use the `CubicSpline` method for demonstration, but the API is the same for all the methods. We can also pass the `extrapolate = true` keyword if we want to allow the interpolation to go beyond the range of the timepoints. The default value is `extrapolate = false`.
 
 ```@example interface
 A1 = CubicSpline(u, t)
@@ -35,6 +35,23 @@ A2(300.0)
     
     The values computed beyond the range of the time points provided during interpolation will not be reliable, as these methods only perform well within the range and the first/last piece polynomial fit is extrapolated on either side which might not reflect the true nature of the data.
 
+The keyword `safetycopy = false` can be passed to make sure no copies of `u` and `t` are made when initializing the interpolation object.
+
+```@example interface
+A3 = QuadraticInterpolation(u, t; safetycopy = false)
+
+# Check for same memory
+u === A3.u.parent
+```
+
+Note that this does not prevent allocation in every interpolation constructor call, because parameter values are cached for all interpolation types except [`ConstantInterpolation`](@ref).
+
+Because of the caching of parameters which depend on `u` and `t`, this data should not be mutated. Therefore `u` and `t` are wrapped in a `ReadOnlyArray` from [ReadOnlyArrays.jl](https://github.com/JuliaArrays/ReadOnlyArrays.jl).
+
+```@repl interface
+A3.t[2] = 3.14
+```
+
 ## Derivatives
 
 Derivatives of the interpolated curves can also be computed at any point for all the methods. Derivatives upto second order is supported where first order derivative is computed analytically and second order using `ForwardDiff.jl`. Order is passed as the third argument. It is 1 by default.
diff --git a/ext/DataInterpolationsOptimExt.jl b/ext/DataInterpolationsOptimExt.jl
index b3bce295..5528503f 100644
--- a/ext/DataInterpolationsOptimExt.jl
+++ b/ext/DataInterpolationsOptimExt.jl
@@ -18,8 +18,9 @@ function Curvefit(u,
         box = false,
         lb = nothing,
         ub = nothing;
-        extrapolate = false)
-    u, t = munge_data(u, t)
+        extrapolate = false,
+        safetycopy = false)
+    u, t = munge_data(u, t, safetycopy)
     errfun(t, u, p) = sum(abs2.(u .- model(t, p)))
     if box == false
         mfit = optimize(p -> errfun(t, u, p), p0, alg)
diff --git a/ext/DataInterpolationsRegularizationToolsExt.jl b/ext/DataInterpolationsRegularizationToolsExt.jl
index a7a914e3..10ea3e4c 100644
--- a/ext/DataInterpolationsRegularizationToolsExt.jl
+++ b/ext/DataInterpolationsRegularizationToolsExt.jl
@@ -69,8 +69,8 @@ A = RegularizationSmooth(u, t, t̂, wls, wr, d; λ = 1.0, alg = :gcv_svd)
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
         wls::AbstractVector, wr::AbstractVector, d::Int = 2;
-        λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     M = _mapping_matrix(t̂, t)
     Wls½ = LA.diagm(sqrt.(wls))
     Wr½ = LA.diagm(sqrt.(wr))
@@ -86,8 +86,8 @@ A = RegularizationSmooth(u, t, d; λ = 1.0, alg = :gcv_svd, extrapolate = false)
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, d::Int = 2;
         λ::Real = 1.0,
-        alg::Symbol = :gcv_svd, extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
@@ -114,8 +114,9 @@ A = RegularizationSmooth(u, t, t̂, d; λ = 1.0, alg = :gcv_svd, extrapolate = f
 ```
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
-        d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd,
+        extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     N, N̂ = length(t), length(t̂)
     M = _mapping_matrix(t̂, t)
     Wls½ = Array{Float64}(LA.I, N, N)
@@ -142,8 +143,8 @@ A = RegularizationSmooth(u, t, t̂, wls, d; λ = 1.0, alg = :gcv_svd, extrapolat
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
         wls::AbstractVector, d::Int = 2; λ::Real = 1.0,
-        alg::Symbol = :gcv_svd, extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     N, N̂ = length(t), length(t̂)
     M = _mapping_matrix(t̂, t)
     Wls½ = LA.diagm(sqrt.(wls))
@@ -171,8 +172,8 @@ A = RegularizationSmooth(
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
         wls::AbstractVector, d::Int = 2; λ::Real = 1.0,
-        alg::Symbol = :gcv_svd, extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
@@ -201,8 +202,8 @@ A = RegularizationSmooth(
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
         wls::AbstractVector, wr::AbstractVector, d::Int = 2;
-        λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
@@ -231,8 +232,8 @@ A = RegularizationSmooth(
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
         wls::Symbol, d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd,
-        extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl
index 7638f1a8..9df9db68 100644
--- a/src/DataInterpolations.jl
+++ b/src/DataInterpolations.jl
@@ -7,6 +7,7 @@ abstract type AbstractInterpolation{T} end
 using LinearAlgebra, RecipesBase
 using PrettyTables
 using ForwardDiff
+using ReadOnlyArrays
 import FindFirstFunctions: searchsortedfirstcorrelated, searchsortedlastcorrelated,
                            bracketstrictlymontonic
 
@@ -74,6 +75,12 @@ function Base.showerror(io::IO, e::DerivativeNotFoundError)
     print(io, DERIVATIVE_NOT_FOUND_ERROR)
 end
 
+const MUST_COPY_ERROR = "A copy must be made of u, t to filter missing data"
+struct MustCopyError <: Exception end
+function Base.showerror(io::IO, e::MustCopyError)
+    print(io, MUST_COPY_ERROR)
+end
+
 export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation,
        AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline,
        BSplineInterpolation, BSplineApprox, CubicHermiteSpline,
@@ -105,12 +112,13 @@ struct RegularizationSmooth{uType, tType, T, T2, ITP <: AbstractInterpolation{T}
             alg,
             Aitp,
             extrapolate)
-        new{typeof(u), typeof(t), eltype(u), typeof(λ), typeof(Aitp)}(u,
-            û,
-            t,
-            t̂,
-            wls,
-            wr,
+        new{typeof(u), typeof(t), eltype(u), typeof(λ), typeof(Aitp)}(
+            readonly_wrap(u),
+            readonly_wrap(û),
+            readonly_wrap(t),
+            readonly_wrap(t̂),
+            readonly_wrap(oftype(u.parent, wls)),
+            readonly_wrap(oftype(u.parent, wr)),
             d,
             λ,
             alg,
diff --git a/src/derivatives.jl b/src/derivatives.jl
index bb3c938e..30c76fd0 100644
--- a/src/derivatives.jl
+++ b/src/derivatives.jl
@@ -17,30 +17,17 @@ function derivative(A, t, order = 1)
     end
 end
 
-function _derivative(A::LinearInterpolation{<:AbstractVector}, t::Number, iguess)
+function _derivative(A::LinearInterpolation, t::Number, iguess)
     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
+    A.p.slope[idx], idx
 end
 
-function _derivative(A::LinearInterpolation{<:AbstractMatrix}, t::Number, iguess)
-    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
-
-function _derivative(A::QuadraticInterpolation{<:AbstractVector}, t::Number, iguess)
-    i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess)
-    dl₀ = (2t - A.t[i₁] - A.t[i₂]) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂]))
-    dl₁ = (2t - A.t[i₀] - A.t[i₂]) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂]))
-    dl₂ = (2t - A.t[i₀] - A.t[i₁]) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁]))
-    A.u[i₀] * dl₀ + A.u[i₁] * dl₁ + A.u[i₂] * dl₂, i₀
-end
-
-function _derivative(A::QuadraticInterpolation{<:AbstractMatrix}, t::Number, iguess)
+function _derivative(A::QuadraticInterpolation, t::Number, iguess)
     i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess)
-    dl₀ = (2t - A.t[i₁] - A.t[i₂]) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂]))
-    dl₁ = (2t - A.t[i₀] - A.t[i₂]) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂]))
-    dl₂ = (2t - A.t[i₀] - A.t[i₁]) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁]))
-    (@views @. A.u[:, i₀] * dl₀ + A.u[:, i₁] * dl₁ + A.u[:, i₂] * dl₂), i₀
+    du₀ = A.p.l₀[i₀] * (2t - A.t[i₁] - A.t[i₂])
+    du₁ = A.p.l₁[i₀] * (2t - A.t[i₀] - A.t[i₂])
+    du₂ = A.p.l₂[i₀] * (2t - A.t[i₀] - A.t[i₁])
+    return @views @. du₀ + du₁ + du₂, i₀
 end
 
 function _derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number)
@@ -125,6 +112,10 @@ function _derivative(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess)
     (@evalpoly wj A.b[idx] 2A.c[j] 3A.d[j]), idx
 end
 
+function _derivative(A::ConstantInterpolation, t::Number, iguess)
+    return zero(first(A.u)), iguess
+end
+
 function _derivative(A::ConstantInterpolation{<:AbstractVector}, t::Number)
     ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError())
     return isempty(searchsorted(A.t, t)) ? zero(A.u[1]) : eltype(A.u)(NaN)
@@ -138,17 +129,18 @@ end
 # QuadraticSpline Interpolation
 function _derivative(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess)
     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.p.σ[idx - 1]
     A.z[idx - 1] + 2σ * (t - A.t[idx - 1]), idx
 end
 
 # CubicSpline Interpolation
 function _derivative(A::CubicSpline{<:AbstractVector}, t::Number, iguess)
     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)
+    Δt₁ = t - A.t[idx]
+    Δt₂ = A.t[idx + 1] - t
+    dI = (-A.z[idx] * Δt₂^2 + A.z[idx + 1] * Δt₁^2) / (2A.h[idx + 1])
+    dC = A.p.c₁[idx]
+    dD = -A.p.c₂[idx]
     dI + dC + dD, idx
 end
 
@@ -160,7 +152,8 @@ function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Num
     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
-    N = DataInterpolations.spline_coefficients(n, A.d - 1, A.k, t_)
+    N = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.N
+    spline_coefficients!(N, A.d - 1, A.k, t_)
     ducum = zero(eltype(A.u))
     if t == A.t[1]
         ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2])
@@ -180,7 +173,8 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig
     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_)
+    N = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.N
+    spline_coefficients!(N, A.d - 1, A.k, t_)
     ducum = zero(eltype(A.u))
     if t == A.t[1]
         ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2])
diff --git a/src/integrals.jl b/src/integrals.jl
index 00036c0a..e64b7568 100644
--- a/src/integrals.jl
+++ b/src/integrals.jl
@@ -25,11 +25,8 @@ end
 function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}},
         idx::Number,
         t::Number)
-    t1 = A.t[idx]
-    t2 = A.t[idx + 1]
-    u1 = A.u[idx]
-    u2 = A.u[idx + 1]
-    t^2 * (u1 - u2) / (2 * t1 - 2 * t2) + t * (t1 * u2 - t2 * u1) / (t1 - t2)
+    Δt = t - A.t[idx]
+    Δt * (A.u[idx] + A.p.slope[idx] * Δt / 2)
 end
 
 function _integral(A::ConstantInterpolation{<:AbstractVector}, idx::Number, t::Number)
@@ -47,49 +44,30 @@ function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}},
         t::Number)
     A.mode == :Backward && idx > 1 && (idx -= 1)
     idx = min(length(A.t) - 2, idx)
-    t1 = A.t[idx]
-    t2 = A.t[idx + 1]
-    t3 = A.t[idx + 2]
-    u1 = A.u[idx]
-    u2 = A.u[idx + 1]
-    u3 = A.u[idx + 2]
-    (t^3 * (-t1 * u2 + t1 * u3 + t2 * u1 - t2 * u3 - t3 * u1 + t3 * u2) /
-     (3 * t1^2 * t2 - 3 * t1^2 * t3 - 3 * t1 * t2^2 + 3 * t1 * t3^2 + 3 * t2^2 * t3 -
-      3 * t2 * t3^2) +
-     t^2 * (t1^2 * u2 - t1^2 * u3 - t2^2 * u1 + t2^2 * u3 + t3^2 * u1 - t3^2 * u2) /
-     (2 * t1^2 * t2 - 2 * t1^2 * t3 - 2 * t1 * t2^2 + 2 * t1 * t3^2 + 2 * t2^2 * t3 -
-      2 * t2 * t3^2) +
-     t *
-     (t1^2 * t2 * u3 - t1^2 * t3 * u2 - t1 * t2^2 * u3 + t1 * t3^2 * u2 + t2^2 * t3 * u1 -
-      t2 * t3^2 * u1) /
-     (t1^2 * t2 - t1^2 * t3 - t1 * t2^2 + t1 * t3^2 + t2^2 * t3 - t2 * t3^2))
+    t₀ = A.t[idx]
+    t₁ = A.t[idx + 1]
+    t₂ = A.t[idx + 2]
+
+    t_sq = (t^2) / 3
+    Iu₀ = A.p.l₀[idx] * t * (t_sq - t * (t₁ + t₂) / 2 + t₁ * t₂)
+    Iu₁ = A.p.l₁[idx] * t * (t_sq - t * (t₀ + t₂) / 2 + t₀ * t₂)
+    Iu₂ = A.p.l₂[idx] * t * (t_sq - t * (t₀ + t₁) / 2 + t₀ * t₁)
+    return Iu₀ + Iu₁ + Iu₂
 end
 
 function _integral(A::QuadraticSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number)
-    t1 = A.t[idx]
-    t2 = A.t[idx + 1]
-    u1 = A.u[idx]
-    z1 = A.z[idx]
-    z2 = A.z[idx + 1]
-    t^3 * (z1 - z2) / (6 * t1 - 6 * t2) + t^2 * (t1 * z2 - t2 * z1) / (2 * t1 - 2 * t2) +
-    t * (-t1^2 * z1 - t1^2 * z2 + 2 * t1 * t2 * z1 + 2 * t1 * u1 - 2 * t2 * u1) /
-    (2 * t1 - 2 * t2)
+    Cᵢ = A.u[idx]
+    Δt = t - A.t[idx]
+    return A.z[idx] * Δt^2 / 2 + A.p.σ[idx] * Δt^3 / 3 + Cᵢ * Δt
 end
 
 function _integral(A::CubicSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number)
-    t1 = A.t[idx]
-    t2 = A.t[idx + 1]
-    u1 = A.u[idx]
-    u2 = A.u[idx + 1]
-    z1 = A.z[idx]
-    z2 = A.z[idx + 1]
-    h2 = A.h[idx + 1]
-    (t^4 * (-z1 + z2) / (24 * h2) + t^3 * (-t1 * z2 + t2 * z1) / (6 * h2) +
-     t^2 * (h2^2 * z1 - h2^2 * z2 + 3 * t1^2 * z2 - 3 * t2^2 * z1 - 6 * u1 + 6 * u2) /
-     (12 * h2) +
-     t *
-     (h2^2 * t1 * z2 - h2^2 * t2 * z1 - t1^3 * z2 - 6 * t1 * u2 + t2^3 * z1 + 6 * t2 * u1) /
-     (6 * h2))
+    Δt₁sq = (t - A.t[idx])^2 / 2
+    Δt₂sq = (A.t[idx + 1] - t)^2 / 2
+    II = (-A.z[idx] * Δt₂sq^2 + A.z[idx + 1] * Δt₁sq^2) / (6A.h[idx + 1])
+    IC = A.p.c₁[idx] * Δt₁sq
+    ID = -A.p.c₂[idx] * Δt₂sq
+    II + IC + ID
 end
 
 function _integral(A::AkimaInterpolation{<:AbstractVector{<:Number}},
@@ -100,24 +78,9 @@ function _integral(A::AkimaInterpolation{<:AbstractVector{<:Number}},
     A.d[idx] * ((t - t1)^4 / 4)
 end
 
-integral(A::LagrangeInterpolation, t1::Number, t2::Number) = throw(IntegralNotFoundError())
-integral(A::LagrangeInterpolation, t::Number) = throw(IntegralNotFoundError())
-
-function integral(A::BSplineInterpolation{<:AbstractVector{<:Number}},
-        t1::Number,
-        t2::Number)
-    throw(IntegralNotFoundError())
-end
-function integral(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number)
-    throw(IntegralNotFoundError())
-end
-
-function integral(A::BSplineApprox{<:AbstractVector{<:Number}}, t1::Number, t2::Number)
-    throw(IntegralNotFoundError())
-end
-function integral(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number)
-    throw(IntegralNotFoundError())
-end
+_integral(A::LagrangeInterpolation, idx::Number, t::Number) = throw(IntegralNotFoundError())
+_integral(A::BSplineInterpolation, idx::Number, t::Number) = throw(IntegralNotFoundError())
+_integral(A::BSplineApprox, idx::Number, t::Number) = throw(IntegralNotFoundError())
 
 # Cubic Hermite Spline
 function _integral(
diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index 15dba20e..bb6339fe 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -12,20 +12,25 @@ Extrapolation extends the last linear polynomial on each side.
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct LinearInterpolation{uType, tType, T} <: AbstractInterpolation{T}
+struct LinearInterpolation{uType, tType, pType, T} <: AbstractInterpolation{T}
     u::uType
     t::tType
+    p::LinearParameterCache{pType}
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
-    function LinearInterpolation(u, t, extrapolate)
-        new{typeof(u), typeof(t), eltype(u)}(u, t, extrapolate, Ref(1))
+    safetycopy::Bool
+    function LinearInterpolation(u, t, p, extrapolate, safetycopy)
+        new{typeof(u), typeof(t), typeof(p.slope), eltype(u)}(
+            u, t, p, extrapolate, Ref(1), safetycopy)
     end
 end
 
-function LinearInterpolation(u, t; extrapolate = false)
-    u, t = munge_data(u, t)
-    LinearInterpolation(u, t, extrapolate)
+function LinearInterpolation(u, t; extrapolate = false, safetycopy = true)
+    u, t = munge_data(u, t, safetycopy)
+    p = LinearParameterCache(u, t)
+    LinearInterpolation(u, t, p, extrapolate, safetycopy)
 end
 
 """
@@ -43,27 +48,34 @@ Extrapolation extends the last quadratic polynomial on each side.
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct QuadraticInterpolation{uType, tType, T} <: AbstractInterpolation{T}
+struct QuadraticInterpolation{uType, tType, pType, T} <: AbstractInterpolation{T}
     u::uType
     t::tType
+    p::QuadraticParameterCache{pType}
     mode::Symbol
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
-    function QuadraticInterpolation(u, t, mode, extrapolate)
+    safetycopy::Bool
+    function QuadraticInterpolation(u, t, p, mode, extrapolate, safetycopy)
         mode ∈ (:Forward, :Backward) ||
             error("mode should be :Forward or :Backward for QuadraticInterpolation")
-        new{typeof(u), typeof(t), eltype(u)}(u, t, mode, extrapolate, Ref(1))
+        new{typeof(u), typeof(t), typeof(p.l₀), eltype(u)}(
+            u, t, p, mode, extrapolate, Ref(1), safetycopy)
     end
 end
 
-function QuadraticInterpolation(u, t, mode; extrapolate = false)
-    u, t = munge_data(u, t)
-    QuadraticInterpolation(u, t, mode, extrapolate)
+function QuadraticInterpolation(u, t, mode; extrapolate = false, safetycopy = true)
+    u, t = munge_data(u, t, safetycopy)
+    p = QuadraticParameterCache(u, t)
+    QuadraticInterpolation(u, t, p, mode, extrapolate, safetycopy)
 end
 
-function QuadraticInterpolation(u, t; extrapolate = false)
-    QuadraticInterpolation(u, t, :Forward; extrapolate)
+function QuadraticInterpolation(u, t; extrapolate = false, safetycopy = true)
+    u, t = munge_data(u, t, safetycopy)
+    p = QuadraticParameterCache(u, t)
+    QuadraticInterpolation(u, t, p, :Forward, extrapolate, safetycopy)
 end
 
 """
@@ -80,6 +92,7 @@ It is the method of interpolation using Lagrange polynomials of (k-1)th order pa
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
 struct LagrangeInterpolation{uType, tType, T, bcacheType} <:
        AbstractInterpolation{T}
@@ -87,27 +100,33 @@ struct LagrangeInterpolation{uType, tType, T, bcacheType} <:
     t::tType
     n::Int
     bcache::bcacheType
+    idxs::Vector{Int}
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
-    function LagrangeInterpolation(u, t, n, extrapolate)
+    safetycopy::Bool
+    function LagrangeInterpolation(u, t, n, extrapolate, safetycopy)
         bcache = zeros(eltype(u[1]), n + 1)
+        idxs = zeros(Int, n + 1)
         fill!(bcache, NaN)
         new{typeof(u), typeof(t), eltype(u), typeof(bcache)}(u,
             t,
             n,
             bcache,
+            idxs,
             extrapolate,
-            Ref(1)
+            Ref(1),
+            safetycopy
         )
     end
 end
 
-function LagrangeInterpolation(u, t, n = length(t) - 1; extrapolate = false)
-    u, t = munge_data(u, t)
+function LagrangeInterpolation(
+        u, t, n = length(t) - 1; extrapolate = false, safetycopy = true)
+    u, t = munge_data(u, t, safetycopy)
     if n != length(t) - 1
         error("Currently only n=length(t) - 1 is supported")
     end
-    LagrangeInterpolation(u, t, n, extrapolate)
+    LagrangeInterpolation(u, t, n, extrapolate, safetycopy)
 end
 
 """
@@ -124,6 +143,7 @@ Extrapolation extends the last cubic polynomial on each side.
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
 struct AkimaInterpolation{uType, tType, bType, cType, dType, T} <:
        AbstractInterpolation{T}
@@ -134,7 +154,8 @@ struct AkimaInterpolation{uType, tType, bType, cType, dType, T} <:
     d::dType
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
-    function AkimaInterpolation(u, t, b, c, d, extrapolate)
+    safetycopy::Bool
+    function AkimaInterpolation(u, t, b, c, d, extrapolate, safetycopy)
         new{typeof(u), typeof(t), typeof(b), typeof(c),
             typeof(d), eltype(u)}(u,
             t,
@@ -142,13 +163,14 @@ struct AkimaInterpolation{uType, tType, bType, cType, dType, T} <:
             c,
             d,
             extrapolate,
-            Ref(1)
+            Ref(1),
+            safetycopy
         )
     end
 end
 
-function AkimaInterpolation(u, t; extrapolate = false)
-    u, t = munge_data(u, t)
+function AkimaInterpolation(u, t; extrapolate = false, safetycopy = true)
+    u, t = munge_data(u, t, safetycopy)
     n = length(t)
     dt = diff(t)
     m = Array{eltype(u)}(undef, n + 3)
@@ -169,7 +191,7 @@ function AkimaInterpolation(u, t; extrapolate = false)
     c = (3.0 .* m[3:(end - 2)] .- 2.0 .* b[1:(end - 1)] .- b[2:end]) ./ dt
     d = (b[1:(end - 1)] .+ b[2:end] .- 2.0 .* m[3:(end - 2)]) ./ dt .^ 2
 
-    AkimaInterpolation(u, t, b, c, d, extrapolate)
+    AkimaInterpolation(u, t, b, c, d, extrapolate, safetycopy)
 end
 
 """
@@ -188,21 +210,25 @@ Extrapolation extends the last constant polynomial at the end points on each sid
 
   - `dir`: indicates which value should be used for interpolation (`:left` or `:right`).
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
 struct ConstantInterpolation{uType, tType, dirType, T} <: AbstractInterpolation{T}
     u::uType
     t::tType
+    p::Nothing
     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, Ref(1))
+    safetycopy::Bool
+    function ConstantInterpolation(u, t, dir, extrapolate, safetycopy)
+        new{typeof(u), typeof(t), typeof(dir), eltype(u)}(
+            u, t, nothing, dir, extrapolate, Ref(1), safetycopy)
     end
 end
 
-function ConstantInterpolation(u, t; dir = :left, extrapolate = false)
-    u, t = munge_data(u, t)
-    ConstantInterpolation(u, t, dir, extrapolate)
+function ConstantInterpolation(u, t; dir = :left, extrapolate = false, safetycopy = true)
+    u, t = munge_data(u, t, safetycopy)
+    ConstantInterpolation(u, t, dir, extrapolate, safetycopy)
 end
 
 """
@@ -219,33 +245,38 @@ Extrapolation extends the last quadratic polynomial on each side.
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct QuadraticSpline{uType, tType, tAType, dType, zType, T} <:
+struct QuadraticSpline{uType, tType, pType, tAType, dType, zType, T} <:
        AbstractInterpolation{T}
     u::uType
     t::tType
+    p::QuadraticSplineParameterCache{pType}
     tA::tAType
     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),
+    safetycopy::Bool
+    function QuadraticSpline(u, t, p, tA, d, z, extrapolate, safetycopy)
+        new{typeof(u), typeof(t), typeof(p.σ), typeof(tA),
             typeof(d), typeof(z), eltype(u)}(u,
             t,
+            p,
             tA,
             d,
             z,
             extrapolate,
-            Ref(1)
+            Ref(1),
+            safetycopy
         )
     end
 end
 
 function QuadraticSpline(u::uType,
         t;
-        extrapolate = false) where {uType <: AbstractVector{<:Number}}
-    u, t = munge_data(u, t)
+        extrapolate = false, safetycopy = true) where {uType <: AbstractVector{<:Number}}
+    u, t = munge_data(u, t, safetycopy)
     s = length(t)
     dl = ones(eltype(t), s - 1)
     d_tmp = ones(eltype(t), s)
@@ -257,11 +288,13 @@ function QuadraticSpline(u::uType,
 
     d = map(i -> i == 1 ? typed_zero : 2 // 1 * (u[i] - u[i - 1]) / (t[i] - t[i - 1]), 1:s)
     z = tA \ d
-    QuadraticSpline(u, t, tA, d, z, extrapolate)
+    p = QuadraticSplineParameterCache(z, t)
+    QuadraticSpline(u, t, p, tA, d, z, extrapolate, safetycopy)
 end
 
-function QuadraticSpline(u::uType, t; extrapolate = false) where {uType <: AbstractVector}
-    u, t = munge_data(u, t)
+function QuadraticSpline(
+        u::uType, t; extrapolate = false, safetycopy = true) where {uType <: AbstractVector}
+    u, t = munge_data(u, t, safetycopy)
     s = length(t)
     dl = ones(eltype(t), s - 1)
     d_tmp = ones(eltype(t), s)
@@ -274,7 +307,8 @@ function QuadraticSpline(u::uType, t; extrapolate = false) where {uType <: Abstr
     d = transpose(reshape(reduce(hcat, d_), :, s))
     z_ = reshape(transpose(tA \ d), size(u[1])..., :)
     z = [z_s for z_s in eachslice(z_, dims = ndims(z_))]
-    QuadraticSpline(u, t, tA, d, z, extrapolate)
+    p = QuadraticSplineParameterCache(z, t)
+    QuadraticSpline(u, t, p, tA, d, z, extrapolate, safetycopy)
 end
 
 """
@@ -291,29 +325,34 @@ Second derivative on both ends are zero, which are also called "natural" boundar
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct CubicSpline{uType, tType, hType, zType, T} <: AbstractInterpolation{T}
+struct CubicSpline{uType, tType, pType, hType, zType, T} <: AbstractInterpolation{T}
     u::uType
     t::tType
+    p::CubicSplineParameterCache{pType}
     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,
+    safetycopy::Bool
+    function CubicSpline(u, t, p, h, z, extrapolate, safetycopy)
+        new{typeof(u), typeof(t), typeof(p.c₁), typeof(h), typeof(z), eltype(u)}(u,
             t,
+            p,
             h,
             z,
             extrapolate,
-            Ref(1)
+            Ref(1),
+            safetycopy
         )
     end
 end
 
 function CubicSpline(u::uType,
         t;
-        extrapolate = false) where {uType <: AbstractVector{<:Number}}
-    u, t = munge_data(u, t)
+        extrapolate = false, safetycopy = true) where {uType <: AbstractVector{<:Number}}
+    u, t = munge_data(u, t, safetycopy)
     n = length(t) - 1
     h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0)
     dl = vcat(h[2:n], zero(eltype(h)))
@@ -330,11 +369,13 @@ function CubicSpline(u::uType,
              6(u[i + 1] - u[i]) / h[i + 1] - 6(u[i] - u[i - 1]) / h[i],
         1:(n + 1))
     z = tA \ d
-    CubicSpline(u, t, h[1:(n + 1)], z, extrapolate)
+    p = CubicSplineParameterCache(u, h, z)
+    CubicSpline(u, t, p, h[1:(n + 1)], z, extrapolate, safetycopy)
 end
 
-function CubicSpline(u::uType, t; extrapolate = false) where {uType <: AbstractVector}
-    u, t = munge_data(u, t)
+function CubicSpline(
+        u::uType, t; extrapolate = false, safetycopy = true) where {uType <: AbstractVector}
+    u, t = munge_data(u, t, safetycopy)
     n = length(t) - 1
     h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0)
     dl = vcat(h[2:n], zero(eltype(h)))
@@ -348,7 +389,8 @@ function CubicSpline(u::uType, t; extrapolate = false) where {uType <: AbstractV
     d = transpose(reshape(reduce(hcat, d_), :, n + 1))
     z_ = reshape(transpose(tA \ d), size(u[1])..., :)
     z = [z_s for z_s in eachslice(z_, dims = ndims(z_))]
-    CubicSpline(u, t, h[1:(n + 1)], z, extrapolate)
+    p = CubicSplineParameterCache(u, h, z)
+    CubicSpline(u, t, p, h[1:(n + 1)], z, extrapolate, safetycopy)
 end
 
 """
@@ -368,8 +410,9 @@ Extrapolation is a constant polynomial of the end points on each side.
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct BSplineInterpolation{uType, tType, pType, kType, cType, T} <:
+struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
        AbstractInterpolation{T}
     u::uType
     t::tType
@@ -377,35 +420,42 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, T} <:
     p::pType  # params vector
     k::kType  # knot vector
     c::cType  # control points
+    N::NType  # Spline coefficients (preallocated memory)
     pVecType::Symbol
     knotVecType::Symbol
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
+    safetycopy::Bool
     function BSplineInterpolation(u,
             t,
             d,
             p,
             k,
             c,
+            N,
             pVecType,
             knotVecType,
-            extrapolate)
-        new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), eltype(u)}(u,
+            extrapolate,
+            safetycopy)
+        new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(N), eltype(u)}(u,
             t,
             d,
             p,
             k,
             c,
+            N,
             pVecType,
             knotVecType,
             extrapolate,
-            Ref(1)
+            Ref(1),
+            safetycopy
         )
     end
 end
 
-function BSplineInterpolation(u, t, d, pVecType, knotVecType; extrapolate = false)
-    u, t = munge_data(u, t)
+function BSplineInterpolation(
+        u, t, d, pVecType, knotVecType; extrapolate = false, safetycopy = true)
+    u, t = munge_data(u, t, safetycopy)
     n = length(t)
     s = zero(eltype(u))
     p = zero(t)
@@ -463,9 +513,12 @@ function BSplineInterpolation(u, t, d, pVecType, knotVecType; extrapolate = fals
         end
     end
     # control points
-    N = spline_coefficients(n, d, k, p)
+    N = zeros(eltype(t), n, n)
+    spline_coefficients!(N, d, k, p)
     c = vec(N \ u[:, :])
-    BSplineInterpolation(u, t, d, p, k, c, pVecType, knotVecType, extrapolate)
+    N = zeros(eltype(t), n)
+    BSplineInterpolation(
+        u, t, d, p, k, c, N, pVecType, knotVecType, extrapolate, safetycopy)
 end
 
 """
@@ -487,8 +540,9 @@ Extrapolation is a constant polynomial of the end points on each side.
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
-struct BSplineApprox{uType, tType, pType, kType, cType, T} <:
+struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
        AbstractInterpolation{T}
     u::uType
     t::tType
@@ -497,10 +551,12 @@ struct BSplineApprox{uType, tType, pType, kType, cType, T} <:
     p::pType  # params vector
     k::kType  # knot vector
     c::cType  # control points
+    N::NType  # Spline coefficients (preallocated memory)
     pVecType::Symbol
     knotVecType::Symbol
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
+    safetycopy::Bool
     function BSplineApprox(u,
             t,
             d,
@@ -508,26 +564,32 @@ struct BSplineApprox{uType, tType, pType, kType, cType, T} <:
             p,
             k,
             c,
+            N,
             pVecType,
             knotVecType,
-            extrapolate)
-        new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), eltype(u)}(u,
+            extrapolate,
+            safetycopy
+    )
+        new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(N), eltype(u)}(u,
             t,
             d,
             h,
             p,
             k,
             c,
+            N,
             pVecType,
             knotVecType,
             extrapolate,
-            Ref(1)
+            Ref(1),
+            safetycopy::Bool
         )
     end
 end
 
-function BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolate = false)
-    u, t = munge_data(u, t)
+function BSplineApprox(
+        u, t, d, h, pVecType, knotVecType; extrapolate = false, safetycopy = true)
+    u, t = munge_data(u, t, safetycopy)
     n = length(t)
     s = zero(eltype(u))
     p = zero(t)
@@ -592,7 +654,7 @@ function BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolate = false)
     q = zeros(eltype(u), n)
     N = zeros(eltype(t), n, h)
     for i in 1:n
-        N[i, :] .= spline_coefficients(h, d, k, p[i])
+        spline_coefficients!(view(N, i, :), d, k, p[i])
     end
     for k in 2:(n - 1)
         q[k] = u[k] - N[k, 1] * u[1] - N[k, h] * u[end]
@@ -609,7 +671,8 @@ function BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolate = false)
     M = transpose(N) * N
     P = M \ Q
     c[2:(end - 1)] .= vec(P)
-    BSplineApprox(u, t, d, h, p, k, c, pVecType, knotVecType, extrapolate)
+    N = zeros(eltype(t), h)
+    BSplineApprox(u, t, d, h, p, k, c, N, pVecType, knotVecType, extrapolate, safetycopy)
 end
 
 """
@@ -626,25 +689,27 @@ It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomi
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
 struct CubicHermiteSpline{uType, tType, duType, pType, T} <: AbstractInterpolation{T}
-    du::uType
+    du::duType
     u::uType
     t::tType
     p::CubicHermiteParameterCache{pType}
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
-    function CubicHermiteSpline(du, u, t, p, extrapolate)
+    safetycopy::Bool
+    function CubicHermiteSpline(du, u, t, p, extrapolate, safetycopy)
         new{typeof(u), typeof(t), typeof(du), typeof(p.c₁), eltype(u)}(
-            du, u, t, p, extrapolate, Ref(1))
+            du, u, t, p, extrapolate, Ref(1), safetycopy)
     end
 end
 
-function CubicHermiteSpline(du, u, t; extrapolate = false)
+function CubicHermiteSpline(du, u, t; extrapolate = false, safetycopy = true)
     @assert length(u)==length(du) "Length of `u` is not equal to length of `du`."
-    u, t = munge_data(u, t)
+    u, t = munge_data(u, t, safetycopy)
     p = CubicHermiteParameterCache(du, u, t)
-    return CubicHermiteSpline(du, u, t, p, extrapolate)
+    return CubicHermiteSpline(du, u, t, p, extrapolate, safetycopy)
 end
 
 """
@@ -662,25 +727,27 @@ It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polyno
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`.
 """
 struct QuinticHermiteSpline{uType, tType, duType, dduType, pType, T} <:
        AbstractInterpolation{T}
-    ddu::uType
-    du::uType
+    ddu::dduType
+    du::duType
     u::uType
     t::tType
     p::QuinticHermiteParameterCache{pType}
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
-    function QuinticHermiteSpline(ddu, du, u, t, p, extrapolate)
+    safetycopy::Bool
+    function QuinticHermiteSpline(ddu, du, u, t, p, extrapolate, safetycopy)
         new{typeof(u), typeof(t), typeof(du), typeof(ddu), typeof(p.c₁), eltype(u)}(
-            ddu, du, u, t, p, extrapolate, Ref(1))
+            ddu, du, u, t, p, extrapolate, Ref(1), safetycopy)
     end
 end
 
-function QuinticHermiteSpline(ddu, du, u, t; extrapolate = false)
+function QuinticHermiteSpline(ddu, du, u, t; extrapolate = false, safetycopy = true)
     @assert length(u)==length(du)==length(ddu) "Length of `u` is not equal to length of `du` or `ddu`."
-    u, t = munge_data(u, t)
+    u, t = munge_data(u, t, safetycopy)
     p = QuinticHermiteParameterCache(ddu, du, u, t)
-    return QuinticHermiteSpline(ddu, du, u, t, p, extrapolate)
+    return QuinticHermiteSpline(ddu, du, u, t, p, extrapolate, safetycopy)
 end
diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl
index 0a543d8f..5d09ceff 100644
--- a/src/interpolation_methods.jl
+++ b/src/interpolation_methods.jl
@@ -13,23 +13,32 @@ function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, igues
         idx = firstindex(A.u) - 1
         t1 = t2 = one(eltype(A.t))
         u1 = u2 = one(eltype(A.u))
+        slope = t * one(eltype(A.p.slope))
     else
         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]
+        slope = A.p.slope[idx]
     end
-    θ = (t - t1) / (t2 - t1)
-    val = (1 - θ) * u1 + θ * u2
-    # Note: The following is limited to when val is NaN as to not change the derivative of exact points.
-    t == t1 && any(isnan, val) && return oftype(val, u1), idx # Return exact value if no interpolation needed (eg when NaN at t2)
-    t == t2 && any(isnan, val) && return oftype(val, u2), idx # ... (eg when NaN at t1)
+
+    Δt = t - t1
+    Δu = slope * Δt
+    val = u1
+    Δu_nan = any(isnan.(Δu))
+    if t == t2 && Δu_nan
+        val = u2
+    elseif !(iszero(Δt) && Δu_nan)
+        val += Δu
+    end
+    val = oftype(Δu, val)
+
     val, idx
 end
 
 function _interpolate(A::LinearInterpolation{<:AbstractMatrix}, t::Number, iguess)
     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
+    Δt = t - A.t[idx]
+    return A.u[:, idx] + A.p.slope[idx] * Δt, idx
 end
 
 # Quadratic Interpolation
@@ -39,87 +48,71 @@ function _quad_interp_indices(A::QuadraticInterpolation, t::Number, iguess)
     idx, idx + 1, idx + 2
 end
 
-function _interpolate(A::QuadraticInterpolation{<:AbstractVector}, t::Number, iguess)
+function _interpolate(A::QuadraticInterpolation, t::Number, iguess)
     i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess)
-    l₀ = ((t - A.t[i₁]) * (t - A.t[i₂])) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂]))
-    l₁ = ((t - A.t[i₀]) * (t - A.t[i₂])) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂]))
-    l₂ = ((t - A.t[i₀]) * (t - A.t[i₁])) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁]))
-    return A.u[i₀] * l₀ + A.u[i₁] * l₁ + A.u[i₂] * l₂, i₀
-end
-
-function _interpolate(A::QuadraticInterpolation{<:AbstractMatrix}, t::Number, iguess)
-    i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess)
-    l₀ = ((t - A.t[i₁]) * (t - A.t[i₂])) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂]))
-    l₁ = ((t - A.t[i₀]) * (t - A.t[i₂])) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂]))
-    l₂ = ((t - A.t[i₀]) * (t - A.t[i₁])) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁]))
-    return A.u[:, i₀] * l₀ + A.u[:, i₁] * l₁ + A.u[:, i₂] * l₂, i₀
+    u₀ = A.p.l₀[i₀] * (t - A.t[i₁]) * (t - A.t[i₂])
+    u₁ = A.p.l₁[i₀] * (t - A.t[i₀]) * (t - A.t[i₂])
+    u₂ = A.p.l₂[i₀] * (t - A.t[i₀]) * (t - A.t[i₁])
+    return u₀ + u₁ + u₂, i₀
 end
 
 # Lagrange Interpolation
-function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number)
-    ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError())
-    idxs = findRequiredIdxs(A, t)
-    if A.t[idxs[1]] == t
-        return A.u[idxs[1]]
+function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess)
+    findRequiredIdxs!(A, t, idx)
+    if A.t[A.idxs[1]] == t
+        return A.u[A.idxs[1]], idx
     end
     N = zero(A.u[1])
     D = zero(A.t[1])
     tmp = N
-    for i in 1:length(idxs)
-        if isnan(A.bcache[idxs[i]])
+    for i in 1:length(A.idxs)
+        if isnan(A.bcache[A.idxs[i]])
             mult = one(A.t[1])
             for j in 1:(i - 1)
-                mult *= (A.t[idxs[i]] - A.t[idxs[j]])
+                mult *= (A.t[A.idxs[i]] - A.t[A.idxs[j]])
             end
-            for j in (i + 1):length(idxs)
-                mult *= (A.t[idxs[i]] - A.t[idxs[j]])
+            for j in (i + 1):length(A.idxs)
+                mult *= (A.t[A.idxs[i]] - A.t[A.idxs[j]])
             end
-            A.bcache[idxs[i]] = mult
+            A.bcache[A.idxs[i]] = mult
         else
-            mult = A.bcache[idxs[i]]
+            mult = A.bcache[A.idxs[i]]
         end
-        tmp = inv((t - A.t[idxs[i]]) * mult)
+        tmp = inv((t - A.t[A.idxs[i]]) * mult)
         D += tmp
-        N += (tmp * A.u[idxs[i]])
+        N += (tmp * A.u[A.idxs[i]])
     end
-    N / D
+    N / D, idx
 end
 
-function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number)
-    ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError())
-    idxs = findRequiredIdxs(A, t)
-    if A.t[idxs[1]] == t
-        return A.u[:, idxs[1]]
+function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess)
+    findRequiredIdxs!(A, t, idx)
+    if A.t[A.idxs[1]] == t
+        return A.u[:, A.idxs[1]], idx
     end
     N = zero(A.u[:, 1])
     D = zero(A.t[1])
     tmp = D
-    for i in 1:length(idxs)
-        if isnan(A.bcache[idxs[i]])
+    for i in 1:length(A.idxs)
+        if isnan(A.bcache[A.idxs[i]])
             mult = one(A.t[1])
             for j in 1:(i - 1)
-                mult *= (A.t[idxs[i]] - A.t[idxs[j]])
+                mult *= (A.t[A.idxs[i]] - A.t[A.idxs[j]])
             end
-            for j in (i + 1):length(idxs)
-                mult *= (A.t[idxs[i]] - A.t[idxs[j]])
+            for j in (i + 1):length(A.idxs)
+                mult *= (A.t[A.idxs[i]] - A.t[A.idxs[j]])
             end
-            A.bcache[idxs[i]] = mult
+            A.bcache[A.idxs[i]] = mult
         else
-            mult = A.bcache[idxs[i]]
+            mult = A.bcache[A.idxs[i]]
         end
-        tmp = inv((t - A.t[idxs[i]]) * mult)
+        tmp = inv((t - A.t[A.idxs[i]]) * mult)
         D += tmp
-        @. N += (tmp * A.u[:, idxs[i]])
+        @. N += (tmp * A.u[:, A.idxs[i]])
     end
-    N / D
-end
-
-function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, idx)
-    _interpolate(A, t), idx
-end
-
-function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, idx)
-    _interpolate(A, t), idx
+    N / D, idx
 end
 
 function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess)
@@ -153,19 +146,20 @@ end
 
 # QuadraticSpline Interpolation
 function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess)
-    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
+    idx = get_idx(A.t, t, iguess)
+    Cᵢ = A.u[idx]
+    Δt = t - A.t[idx]
+    return A.z[idx] * Δt + A.p.σ[idx] * Δt^2 + Cᵢ, idx
 end
 
 # CubicSpline Interpolation
 function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess)
     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)
+    Δt₁ = t - A.t[idx]
+    Δt₂ = A.t[idx + 1] - t
+    I = (A.z[idx] * Δt₂^3 + A.z[idx + 1] * Δt₁^3) / (6A.h[idx + 1])
+    C = A.p.c₁[idx] * Δt₁
+    D = A.p.c₂[idx] * Δt₂
     I + C + D, idx
 end
 
@@ -179,9 +173,10 @@ function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}},
     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)
+    N = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.N
+    nonzero_coefficient_idxs = spline_coefficients!(N, A.d, A.k, t)
     ucum = zero(eltype(A.u))
-    for i in 1:n
+    for i in nonzero_coefficient_idxs
         ucum += N[i] * A.c[i]
     end
     ucum, idx
@@ -194,9 +189,10 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i
     # change t into param [0 1]
     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 = spline_coefficients(A.h, A.d, A.k, t)
+    N = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.N
+    nonzero_coefficient_idxs = spline_coefficients!(N, A.d, A.k, t)
     ucum = zero(eltype(A.u))
-    for i in 1:(A.h)
+    for i in nonzero_coefficient_idxs
         ucum += N[i] * A.c[i]
     end
     ucum, idx
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index bb62e78a..370c4062 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -1,14 +1,42 @@
-function findRequiredIdxs(A::LagrangeInterpolation, t)
-    idxs = sortperm(A.t, by = x -> abs(t - x))
-    idxs[1:(A.n + 1)]
+function findRequiredIdxs!(A::LagrangeInterpolation, t, idx)
+    n = length(A.t) - 1
+    i_min, idx_min, idx_max = if t == A.t[idx]
+        A.idxs[1] = idx
+        2, idx, idx
+    else
+        1, idx + 1, idx
+    end
+    for i in i_min:(n + 1)
+        if idx_min == 1
+            A.idxs[i:end] .= range(idx_max + 1, idx_max + (n + 2 - i))
+            break
+        elseif idx_max == length(A.t)
+            A.idxs[i:end] .= (idx_min - 1):-1:(idx_min - (n + 2 - i))
+            break
+        else
+            left_diff = abs(t - A.t[idx_min - 1])
+            right_diff = abs(t - A.t[idx_max + 1])
+            left_expand = left_diff <= right_diff
+        end
+        if left_expand
+            idx_min -= 1
+            A.idxs[i] = idx_min
+        else
+            idx_max += 1
+            A.idxs[i] = idx_max
+        end
+    end
+    return idx
 end
 
-function spline_coefficients(n, d, k, u::Number)
-    N = zeros(eltype(u), n)
+function spline_coefficients!(N, d, k, u::Number)
+    N .= zero(u)
     if u == k[1]
         N[1] = one(u)
+        return 1:1
     elseif u == k[end]
         N[end] = one(u)
+        return length(N):length(N)
     else
         i = findfirst(x -> x > u, k) - 1
         N[i] = one(u)
@@ -20,51 +48,68 @@ function spline_coefficients(n, d, k, u::Number)
             end
             N[i] = (u - k[i]) / (k[i + deg] - k[i]) * N[i]
         end
+        return (i - d):i
     end
-    N
 end
 
-function spline_coefficients(n, d, k, u::AbstractVector)
-    N = zeros(eltype(u), n, n)
-    for i in 1:n
-        N[i, :] .= spline_coefficients(n, d, k, u[i])
+function spline_coefficients!(N, d, k, u::AbstractVector)
+    for i in 1:size(N)[2]
+        spline_coefficients!(view(N, i, :), d, k, u[i])
     end
-    N
+    return nothing
 end
 
 # helper function for data manipulation
-function munge_data(u::AbstractVector{<:Real}, t::AbstractVector{<:Real})
-    return u, t
+function munge_data(u::AbstractVector{<:Real}, t::AbstractVector{<:Real}, safetycopy::Bool)
+    if safetycopy
+        u = copy(u)
+        t = copy(t)
+    end
+    return readonly_wrap(u), readonly_wrap(t)
 end
 
-function munge_data(u::AbstractVector, t::AbstractVector)
+function munge_data(u::AbstractVector, t::AbstractVector, safetycopy::Bool)
     Tu = Base.nonmissingtype(eltype(u))
     Tt = Base.nonmissingtype(eltype(t))
     @assert length(t) == length(u)
     non_missing_indices = collect(
-        idx for idx in 1:length(t)
-    if !ismissing(u[idx]) && !ismissing(t[idx])
+        i for i in 1:length(t)
+    if !ismissing(u[i]) && !ismissing(t[i])
     )
-    newu = Tu.([u[idx] for idx in non_missing_indices])
-    newt = Tt.([t[idx] for idx in non_missing_indices])
 
-    return newu, newt
+    if safetycopy
+        u = Tu.([u[i] for i in non_missing_indices])
+        t = Tt.([t[i] for i in non_missing_indices])
+    else
+        !isempty(non_missing_indices) && throw(MustCopyError())
+    end
+
+    return readonly_wrap(u), readonly_wrap(t)
 end
 
-function munge_data(U::StridedMatrix, t::AbstractVector)
+function munge_data(U::StridedMatrix, t::AbstractVector, safetycopy::Bool)
     TU = Base.nonmissingtype(eltype(U))
     Tt = Base.nonmissingtype(eltype(t))
     @assert length(t) == size(U, 2)
     non_missing_indices = collect(
-        idx for idx in 1:length(t)
-    if !any(ismissing, U[:, idx]) && !ismissing(t[idx])
+        i for i in 1:length(t)
+    if !any(ismissing, U[:, i]) && !ismissing(t[i])
     )
-    newUs = [TU.(U[:, idx]) for idx in non_missing_indices]
-    newt = Tt.([t[idx] for idx in non_missing_indices])
 
-    return hcat(newUs...), newt
+    if safetycopy
+        U = hcat([TU.(U[:, i]) for i in non_missing_indices]...)
+        t = Tt.([t[i] for i in non_missing_indices])
+    else
+        !isempty(non_missing_indices) && throw(MustCopyError())
+    end
+
+    return readonly_wrap(U), readonly_wrap(t)
 end
 
+# Don't nest ReadOnlyArrays
+readonly_wrap(a::AbstractArray) = ReadOnlyArray(a)
+readonly_wrap(a::ReadOnlyArray) = a
+
 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
diff --git a/src/online.jl b/src/online.jl
index 4ad06f3a..dc500611 100644
--- a/src/online.jl
+++ b/src/online.jl
@@ -1,40 +1,63 @@
 import Base: append!, push!
 
-function push!(di::LinearInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T}
-    push!(di.u, u)
-    push!(di.t, t)
-    di
+function push!(A::LinearInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T}
+    push!(A.u.parent, u)
+    push!(A.t.parent, t)
+    slope = linear_interpolation_parameters(A.u, A.t, length(A.t) - 1)
+    push!(A.p.slope, slope)
+    A
 end
 
-function push!(di::QuadraticInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T}
-    push!(di.u, u)
-    push!(di.t, t)
-    di
+function push!(A::QuadraticInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T}
+    push!(A.u.parent, u)
+    push!(A.t.parent, t)
+    l₀, l₁, l₂ = quadratic_interpolation_parameters(A.u, A.t, length(A.t) - 2)
+    push!(A.p.l₀, l₀)
+    push!(A.p.l₁, l₁)
+    push!(A.p.l₂, l₂)
+    A
 end
 
-function push!(di::ConstantInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T}
-    push!(di.u, u)
-    push!(di.t, t)
-    di
+function push!(A::ConstantInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T}
+    push!(A.u.parent, u)
+    push!(A.t.parent, t)
+    A
 end
 
-function append!(di::LinearInterpolation{U, T}, u::U, t::T) where {U, T}
-    u, t = munge_data(u, t)
-    append!(di.u, u)
-    append!(di.t, t)
-    di
+function append!(
+        A::LinearInterpolation{ReadOnlyVector{eltypeU, U}, ReadOnlyVector{eltypeT, T}}, u::U, t::T) where {
+        eltypeU, U, eltypeT, T}
+    length_old = length(A.t)
+    u, t = munge_data(u, t, true)
+    append!(A.u.parent, u)
+    append!(A.t.parent, t)
+    slope = linear_interpolation_parameters.(
+        Ref(A.u), Ref(A.t), length_old:(length(A.t) - 1))
+    append!(A.p.slope, slope)
+    A
 end
 
-function append!(di::QuadraticInterpolation{U, T}, u::U, t::T) where {U, T}
-    u, t = munge_data(u, t)
-    append!(di.u, u)
-    append!(di.t, t)
-    di
+function append!(
+        A::ConstantInterpolation{ReadOnlyVector{eltypeU, U}, ReadOnlyVector{eltypeT, T}}, u::U, t::T) where {
+        eltypeU, U, eltypeT, T}
+    u, t = munge_data(u, t, true)
+    append!(A.u.parent, u)
+    append!(A.t.parent, t)
+    A
 end
 
-function append!(di::ConstantInterpolation{U, T}, u::U, t::T) where {U, T}
-    u, t = munge_data(u, t)
-    append!(di.u, u)
-    append!(di.t, t)
-    di
+function append!(
+        A::QuadraticInterpolation{ReadOnlyVector{eltypeU, U}, ReadOnlyVector{eltypeT, T}}, u::U, t::T) where {
+        eltypeU, U, eltypeT, T}
+    length_old = length(A.t)
+    u, t = munge_data(u, t, true)
+    append!(A.u.parent, u)
+    append!(A.t.parent, t)
+    parameters = quadratic_interpolation_parameters.(
+        Ref(A.u), Ref(A.t), (length_old - 1):(length(A.t) - 2))
+    l₀, l₁, l₂ = collect.(eachrow(hcat(collect.(parameters)...)))
+    append!(A.p.l₀, l₀)
+    append!(A.p.l₁, l₁)
+    append!(A.p.l₂, l₂)
+    A
 end
diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl
index 66aa218f..2820dc8f 100644
--- a/src/parameter_caches.jl
+++ b/src/parameter_caches.jl
@@ -1,16 +1,100 @@
+struct LinearParameterCache{pType}
+    slope::pType
+end
+
+function LinearParameterCache(u, t)
+    slope = linear_interpolation_parameters.(Ref(u), Ref(t), 1:(length(t) - 1))
+    return LinearParameterCache(slope)
+end
+
+function linear_interpolation_parameters(u, t, idx)
+    Δu = u isa AbstractMatrix ? u[:, idx + 1] - u[:, idx] : u[idx + 1] - u[idx]
+    Δt = t[idx + 1] - t[idx]
+    slope = Δu / Δt
+    slope = iszero(Δt) ? zero(slope) : slope
+    return slope
+end
+
+struct QuadraticParameterCache{pType}
+    l₀::pType
+    l₁::pType
+    l₂::pType
+end
+
+function QuadraticParameterCache(u, t)
+    parameters = quadratic_interpolation_parameters.(
+        Ref(u), Ref(t), 1:(length(t) - 2))
+    l₀, l₁, l₂ = collect.(eachrow(hcat(collect.(parameters)...)))
+    return QuadraticParameterCache(l₀, l₁, l₂)
+end
+
+function quadratic_interpolation_parameters(u, t, idx)
+    if u isa AbstractMatrix
+        u₀ = u[:, idx]
+        u₁ = u[:, idx + 1]
+        u₂ = u[:, idx + 2]
+    else
+        u₀ = u[idx]
+        u₁ = u[idx + 1]
+        u₂ = u[idx + 2]
+    end
+    t₀ = t[idx]
+    t₁ = t[idx + 1]
+    t₂ = t[idx + 2]
+    Δt₀ = t₁ - t₀
+    Δt₁ = t₂ - t₁
+    Δt₂ = t₂ - t₀
+    l₀ = u₀ / (Δt₀ * Δt₂)
+    l₁ = -u₁ / (Δt₀ * Δt₁)
+    l₂ = u₂ / (Δt₂ * Δt₁)
+    return l₀, l₁, l₂
+end
+
+struct QuadraticSplineParameterCache{pType}
+    σ::pType
+end
+
+function QuadraticSplineParameterCache(z, t)
+    σ = quadratic_spline_parameters.(Ref(z), Ref(t), 1:(length(t) - 1))
+    return QuadraticSplineParameterCache(σ)
+end
+
+function quadratic_spline_parameters(z, t, idx)
+    σ = 1 // 2 * (z[idx + 1] - z[idx]) / (t[idx + 1] - t[idx])
+    return σ
+end
+
+struct CubicSplineParameterCache{pType}
+    c₁::pType
+    c₂::pType
+end
+
+function CubicSplineParameterCache(u, h, z)
+    parameters = cubic_spline_parameters.(
+        Ref(u), Ref(h), Ref(z), 1:(size(u)[end] - 1))
+    c₁, c₂ = collect.(eachrow(hcat(collect.(parameters)...)))
+    return CubicSplineParameterCache(c₁, c₂)
+end
+
+function cubic_spline_parameters(u, h, z, idx)
+    c₁ = (u[idx + 1] / h[idx + 1] - z[idx + 1] * h[idx + 1] / 6)
+    c₂ = (u[idx] / h[idx + 1] - z[idx] * h[idx + 1] / 6)
+    return c₁, c₂
+end
+
 struct CubicHermiteParameterCache{pType}
     c₁::pType
     c₂::pType
 end
 
 function CubicHermiteParameterCache(du, u, t)
-    parameters = CubicHermiteSplineParameters.(
+    parameters = cubic_hermite_spline_parameters.(
         Ref(du), Ref(u), Ref(t), 1:(length(t) - 1))
     c₁, c₂ = collect.(eachrow(hcat(collect.(parameters)...)))
     return CubicHermiteParameterCache(c₁, c₂)
 end
 
-function CubicHermiteSplineParameters(du, u, t, idx)
+function cubic_hermite_spline_parameters(du, u, t, idx)
     Δt = t[idx + 1] - t[idx]
     u₀ = u[idx]
     u₁ = u[idx + 1]
@@ -28,13 +112,13 @@ struct QuinticHermiteParameterCache{pType}
 end
 
 function QuinticHermiteParameterCache(ddu, du, u, t)
-    parameters = QuinticHermiteSplineParameters.(
+    parameters = quintic_hermite_spline_parameters.(
         Ref(ddu), Ref(du), Ref(u), Ref(t), 1:(length(t) - 1))
     c₁, c₂, c₃ = collect.(eachrow(hcat(collect.(parameters)...)))
     return QuinticHermiteParameterCache(c₁, c₂, c₃)
 end
 
-function QuinticHermiteSplineParameters(ddu, du, u, t, idx)
+function quintic_hermite_spline_parameters(ddu, du, u, t, idx)
     Δt = t[idx + 1] - t[idx]
     u₀ = u[idx]
     u₁ = u[idx + 1]
diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl
index 225ad9f3..37351d0d 100644
--- a/test/derivative_tests.jl
+++ b/test/derivative_tests.jl
@@ -125,6 +125,13 @@ end
     end
 end
 
+@testset "Constant Interpolation" begin
+    u = [0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0]
+    t = collect(0.0:10.0)
+    A = ConstantInterpolation(u, t)
+    @test all(derivative.(Ref(A), t) .== 0.0)
+end
+
 @testset "Quadratic Spline" begin
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
diff --git a/test/integral_tests.jl b/test/integral_tests.jl
index a5641dff..3d59a6c0 100644
--- a/test/integral_tests.jl
+++ b/test/integral_tests.jl
@@ -83,7 +83,7 @@ end
 
 @testset "LagrangeInterpolation" begin
     u = [1.0, 4.0, 9.0]
-    t = [1.0, 2.0, 3.0]
+    t = [1.0, 2.0, 6.0]
     A = LagrangeInterpolation(u, t)
     @test_throws DataInterpolations.IntegralNotFoundError integral(A, 1.0, 2.0)
     @test_throws DataInterpolations.IntegralNotFoundError integral(A, 5.0)
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index f9ef87f8..779c5bd8 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -3,6 +3,18 @@ using FindFirstFunctions: searchsortedfirstcorrelated
 using StableRNGs
 using Optim, ForwardDiff
 
+function test_interpolation_type(T)
+    @test T <: DataInterpolations.AbstractInterpolation
+    @test hasfield(T, :u)
+    @test hasfield(T, :t)
+    @test hasfield(T, :extrapolate)
+    @test hasfield(T, :idx_prev)
+    @test hasfield(T, :safetycopy)
+    @test !isempty(methods(DataInterpolations._interpolate, (T, Any, Number)))
+    @test !isempty(methods(DataInterpolations._integral, (T, Any, Number)))
+    @test !isempty(methods(DataInterpolations._derivative, (T, Any, Number)))
+end
+
 function test_cached_index(A)
     for t in range(first(A.t), last(A.t); length = 2 * length(A.t) - 1)
         A(t)
@@ -13,6 +25,8 @@ function test_cached_index(A)
 end
 
 @testset "Linear Interpolation" begin
+    test_interpolation_type(LinearInterpolation)
+
     for t in (1.0:10.0, 1.0collect(1:10))
         u = 2.0collect(1:10)
         #t = 1.0collect(1:10)
@@ -75,14 +89,6 @@ end
     @test A(3.0) == 2.0
     @test A(4.0) == 3.0
 
-    u = [0.0, NaN, 2.0, 3.0]
-    A = LinearInterpolation(u, t; extrapolate = true)
-    @test A(1.0) == 0.0
-    @test isnan(A(2.0))
-    @test isnan(A(2.5))
-    @test A(3.0) == 2.0
-    @test A(4.0) == 3.0
-
     u = [0.0, 1.0, NaN, 3.0]
     A = LinearInterpolation(u, t; extrapolate = true)
     @test A(1.0) == 0.0
@@ -167,6 +173,8 @@ end
 end
 
 @testset "Quadratic Interpolation" begin
+    test_interpolation_type(QuadraticInterpolation)
+
     u = [1.0, 4.0, 9.0, 16.0]
     t = [1.0, 2.0, 3.0, 4.0]
     A = QuadraticInterpolation(u, t; extrapolate = true)
@@ -262,6 +270,8 @@ end
 end
 
 @testset "Lagrange Interpolation" begin
+    test_interpolation_type(LagrangeInterpolation)
+
     u = [1.0, 4.0, 9.0]
     t = [1.0, 2.0, 3.0]
     A = LagrangeInterpolation(u, t)
@@ -320,6 +330,8 @@ end
 end
 
 @testset "Akima Interpolation" begin
+    test_interpolation_type(AkimaInterpolation)
+
     u = [0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0]
     t = collect(0.0:10.0)
     A = AkimaInterpolation(u, t)
@@ -349,6 +361,8 @@ end
 end
 
 @testset "ConstantInterpolation" begin
+    test_interpolation_type(ConstantInterpolation)
+
     t = [1.0, 2.0, 3.0, 4.0]
 
     @testset "Vector case" for u in [[1.0, 2.0, 0.0, 1.0], ["B", "C", "A", "B"]]
@@ -473,6 +487,8 @@ end
 end
 
 @testset "QuadraticSpline Interpolation" begin
+    test_interpolation_type(QuadraticSpline)
+
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
 
@@ -518,6 +534,8 @@ end
 end
 
 @testset "CubicSpline Interpolation" begin
+    test_interpolation_type(CubicSpline)
+
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
 
@@ -569,11 +587,13 @@ end
 end
 
 @testset "BSplines" begin
+
     # BSpline Interpolation and Approximation
     t = [0, 62.25, 109.66, 162.66, 205.8, 252.3]
     u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
 
     @testset "BSplineInterpolation" begin
+        test_interpolation_type(BSplineInterpolation)
         A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform)
 
         @test [A(25.0), A(80.0)] == [13.454197730061425, 10.305633616059845]
@@ -605,6 +625,7 @@ end
     end
 
     @testset "BSplineApprox" begin
+        test_interpolation_type(BSplineApprox)
         A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform)
 
         @test [A(25.0), A(80.0)] ≈ [12.979802931218234, 10.914310609953178]
@@ -623,6 +644,8 @@ end
 end
 
 @testset "Cubic Hermite Spline" begin
+    test_interpolation_type(CubicHermiteSpline)
+
     du = [-0.047, -0.058, 0.054, 0.012, -0.068, 0.0]
     u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
     t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3]
@@ -636,6 +659,8 @@ end
 end
 
 @testset "Quintic Hermite Spline" begin
+    test_interpolation_type(QuinticHermiteSpline)
+
     ddu = [0.0, -0.00033, 0.0051, -0.0067, 0.0029, 0.0]
     du = [-0.047, -0.058, 0.054, 0.012, -0.068, 0.0]
     u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
diff --git a/test/online_tests.jl b/test/online_tests.jl
index cc4d4c84..571c2dc0 100644
--- a/test/online_tests.jl
+++ b/test/online_tests.jl
@@ -1,18 +1,31 @@
 using DataInterpolations
 
-t = [1, 2, 3]
-u = [0, 1, 0]
+t1 = [1.0, 2.0, 3.0]
+u1 = [0.0, 1.0, 0.0]
 
-for di in [LinearInterpolation, QuadraticInterpolation, ConstantInterpolation]
-    li = di(copy(u), copy(t))
-    append!(li, u, t)
-    li2 = di(vcat(u, u), vcat(t, t))
-    @test li.u == li2.u
-    @test li.t == li2.t
+t2 = [4.0, 5.0, 6.0]
+u2 = [1.0, 2.0, 1.0]
 
-    li = di(copy(u), copy(t))
-    push!(li, 1, 4)
-    li2 = di(vcat(u, 1), vcat(t, 4))
-    @test li.u == li2.u
-    @test li.t == li2.t
+ts = 1.0:0.5:6.0
+
+for method in [LinearInterpolation, QuadraticInterpolation, ConstantInterpolation]
+    func1 = method(u1, t1)
+    append!(func1, u2, t2)
+    func2 = method(vcat(u1, u2), vcat(t1, t2))
+    @test func1.u == func2.u
+    @test func1.t == func2.t
+    for name in propertynames(func1.p)
+        @test getfield(func1.p, name) == getfield(func2.p, name)
+    end
+    @test func1(ts) == func2(ts)
+
+    func1 = method(u1, t1)
+    push!(func1, 1.0, 4.0)
+    func2 = method(vcat(u1, 1.0), vcat(t1, 4.0))
+    @test func1.u == func2.u
+    @test func1.t == func2.t
+    for name in propertynames(func1.p)
+        @test getfield(func1.p, name) == getfield(func2.p, name)
+    end
+    @test func1(ts) == func2(ts)
 end
diff --git a/test/parameter_tests.jl b/test/parameter_tests.jl
index 3bedf3bc..bcd26cf7 100644
--- a/test/parameter_tests.jl
+++ b/test/parameter_tests.jl
@@ -1,5 +1,36 @@
 using DataInterpolations
 
+@testset "Linear Interpolation" begin
+    u = [1.0, 5.0, 3.0, 4.0, 4.0]
+    t = collect(1:5)
+    A = LinearInterpolation(u, t)
+    @test A.p.slope ≈ [4.0, -2.0, 1.0, 0.0]
+end
+
+@testset "Quadratic Interpolation" begin
+    u = [1.0, 5.0, 3.0, 4.0, 4.0]
+    t = collect(1:5)
+    A = QuadraticInterpolation(u, t)
+    @test A.p.l₀ ≈ [0.5, 2.5, 1.5]
+    @test A.p.l₁ ≈ [-5.0, -3.0, -4.0]
+    @test A.p.l₂ ≈ [1.5, 2.0, 2.0]
+end
+
+@testset "Quadratic Spline" begin
+    u = [1.0, 5.0, 3.0, 4.0, 4.0]
+    t = collect(1:5)
+    A = QuadraticSpline(u, t)
+    @test A.p.σ ≈ [4.0, -10.0, 13.0, -14.0]
+end
+
+@testset "Cubic Spline" begin
+    u = [1, 5, 3, 4, 4]
+    t = collect(1:5)
+    A = CubicSpline(u, t)
+    @test A.p.c₁ ≈ [6.839285714285714, 1.642857142857143, 4.589285714285714, 4.0]
+    @test A.p.c₂ ≈ [1.0, 6.839285714285714, 1.642857142857143, 4.589285714285714]
+end
+
 @testset "Cubic Hermite Spline" begin
     du = [5.0, 3.0, 6.0, 8.0, 1.0]
     u = [1.0, 5.0, 3.0, 4.0, 4.0]