diff --git a/README.md b/README.md
index c3ac44c8..8a326579 100644
--- a/README.md
+++ b/README.md
@@ -59,6 +59,8 @@ corresponding to `(u,t)` pairs.
       + `pVec` - Symbol to Parameters Vector, `pVec = :Uniform` for uniform spaced parameters and `pVec = :ArcLen` for parameters generated by chord length method.
       + `knotVec` - Symbol to Knot Vector, `knotVec = :Uniform` for uniform knot vector, `knotVec = :Average` for average spaced knot vector.
   - `BSplineApprox(u,t,d,h,pVec,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h<length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
+  - `CubicHermiteSpline(du, u, t)` - A third order Hermite interpolation, which matches the values and first (`du`) order derivatives in the data points exactly.
+  - `QuinticHermiteSpline(ddu, du, u, t)` - A fifth order Hermite interpolation, which matches the values and first (`du`) and second (`ddu`) order derivatives in the data points exactly.
 
 ## Extension Methods
 
diff --git a/docs/src/index.md b/docs/src/index.md
index 7e8d1ee7..f5b1512b 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -34,6 +34,8 @@ corresponding to `(u,t)` pairs.
       + `pVec` - Symbol to Parameters Vector, `pVec = :Uniform` for uniformly spaced parameters, and `pVec = :ArcLen` for parameters generated by the chord length method.
       + `knotVec` - Symbol to Knot Vector, `knotVec = :Uniform` for uniform knot vector, `knotVec = :Average` for average spaced knot vector.
   - `BSplineApprox(u,t,d,h,pVec,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h<length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
+  - `CubicHermiteSpline(du, u, t)` - A third order Hermite interpolation, which matches the values and first (`du`) order derivatives in the data points exactly.
+  - `QuinticHermiteSpline(ddu, du, u, t)` - a fifth order Hermite interpolation, which matches the values and first (`du`) and second (`ddu`) order derivatives in the data points exactly.
 
 ## Extension Methods
 
diff --git a/docs/src/manual.md b/docs/src/manual.md
index 5081b71c..7dddd3af 100644
--- a/docs/src/manual.md
+++ b/docs/src/manual.md
@@ -10,4 +10,6 @@ QuadraticSpline
 CubicSpline
 BSplineInterpolation
 BSplineApprox
+CubicHermiteSpline
+QuinticHermiteSpline
 ```
diff --git a/docs/src/methods.md b/docs/src/methods.md
index 7a8f2ddf..624407f7 100644
--- a/docs/src/methods.md
+++ b/docs/src/methods.md
@@ -134,6 +134,29 @@ scatter(t, u, label = "input data")
 plot!(A)
 ```
 
+## Cubic Hermite Spline
+
+This is the cubic (third order) Hermite interpolation. It matches the values and first order derivatives in the data points exactly.
+
+```@example tutorial
+du = [-0.047, -0.058, 0.054, 0.012, -0.068, 0.0011]
+A = CubicHermiteSpline(du, u, t)
+scatter(t, u, label = "input data")
+plot!(A)
+```
+
+## Quintic Hermite Spline
+
+This is the quintic (fifth order) Hermite interpolation. It matches the values and first and second order derivatives in the data points exactly.
+
+```@example tutorial
+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.0011]
+A = QuinticHermiteSpline(ddu, du, u, t)
+scatter(t, u, label = "input data")
+plot!(A)
+```
+
 ## Regularization Smoothing
 
 Smoothing by regularization (a.k.a. ridge regression) finds a function ``\hat{u}``
diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl
index 38cb387d..7638f1a8 100644
--- a/src/DataInterpolations.jl
+++ b/src/DataInterpolations.jl
@@ -10,6 +10,7 @@ using ForwardDiff
 import FindFirstFunctions: searchsortedfirstcorrelated, searchsortedlastcorrelated,
                            bracketstrictlymontonic
 
+include("parameter_caches.jl")
 include("interpolation_caches.jl")
 include("interpolation_utils.jl")
 include("interpolation_methods.jl")
@@ -75,7 +76,8 @@ end
 
 export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation,
        AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline,
-       BSplineInterpolation, BSplineApprox
+       BSplineInterpolation, BSplineApprox, CubicHermiteSpline,
+       QuinticHermiteSpline
 
 # added for RegularizationSmooth, JJS 11/27/21
 ### Regularization data smoothing and interpolation
diff --git a/src/derivatives.jl b/src/derivatives.jl
index 3fb10356..bb3c938e 100644
--- a/src/derivatives.jl
+++ b/src/derivatives.jl
@@ -191,3 +191,26 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig
     end
     ducum * A.d * scale, idx
 end
+
+# Cubic Hermite Spline
+function _derivative(
+        A::CubicHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess)
+    Δt₀ = t - A.t[idx]
+    Δt₁ = t - A.t[idx + 1]
+    out = A.du[idx]
+    out += Δt₀ * (Δt₀ * A.p.c₂[idx] + 2(A.p.c₁[idx] + Δt₁ * A.p.c₂[idx]))
+    out, idx
+end
+
+# Quintic Hermite Spline
+function _derivative(
+        A::QuinticHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess)
+    Δt₀ = t - A.t[idx]
+    Δt₁ = t - A.t[idx + 1]
+    out = A.du[idx] + A.ddu[idx] * Δt₀
+    out += Δt₀^2 *
+           (3A.p.c₁[idx] + (3Δt₁ + Δt₀) * A.p.c₂[idx] + (3Δt₁^2 + Δt₀ * 2Δt₁) * A.p.c₃[idx])
+    out, idx
+end
diff --git a/src/integrals.jl b/src/integrals.jl
index 2751eb56..00036c0a 100644
--- a/src/integrals.jl
+++ b/src/integrals.jl
@@ -118,3 +118,28 @@ end
 function integral(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number)
     throw(IntegralNotFoundError())
 end
+
+# Cubic Hermite Spline
+function _integral(
+        A::CubicHermiteSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number)
+    Δt₀ = t - A.t[idx]
+    Δt₁ = t - A.t[idx + 1]
+    out = Δt₀ * (A.u[idx] + Δt₀ * A.du[idx] / 2)
+    p = A.p.c₁[idx] + Δt₁ * A.p.c₂[idx]
+    dp = A.p.c₂[idx]
+    out += Δt₀^3 / 3 * (p - dp * Δt₀ / 4)
+    out
+end
+
+# Quintic Hermite Spline
+function _integral(
+        A::QuinticHermiteSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number)
+    Δt₀ = t - A.t[idx]
+    Δt₁ = t - A.t[idx + 1]
+    out = Δt₀ * (A.u[idx] + A.du[idx] * Δt₀ / 2 + A.ddu[idx] * Δt₀^2 / 6)
+    p = A.p.c₁[idx] + A.p.c₂[idx] * Δt₁ + A.p.c₃[idx] * Δt₁^2
+    dp = A.p.c₂[idx] + 2A.p.c₃[idx] * Δt₁
+    ddp = 2A.p.c₃[idx]
+    out += Δt₀^4 / 4 * (p - Δt₀ / 5 * dp + Δt₀^2 / 30 * ddp)
+    out
+end
diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index 86c6be61..15dba20e 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -611,3 +611,76 @@ function BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolate = false)
     c[2:(end - 1)] .= vec(P)
     BSplineApprox(u, t, d, h, p, k, c, pVecType, knotVecType, extrapolate)
 end
+
+"""
+    CubicHermiteSpline(du, u, t; extrapolate = false)
+
+It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomial such that the value and the first derivative are equal to given values in the data points.
+
+## Arguments
+
+  - `du`: the derivative at the data points.
+  - `u`: data points.
+  - `t`: time points.
+
+## Keyword Arguments
+
+  - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+"""
+struct CubicHermiteSpline{uType, tType, duType, pType, T} <: AbstractInterpolation{T}
+    du::uType
+    u::uType
+    t::tType
+    p::CubicHermiteParameterCache{pType}
+    extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
+    function CubicHermiteSpline(du, u, t, p, extrapolate)
+        new{typeof(u), typeof(t), typeof(du), typeof(p.c₁), eltype(u)}(
+            du, u, t, p, extrapolate, Ref(1))
+    end
+end
+
+function CubicHermiteSpline(du, u, t; extrapolate = false)
+    @assert length(u)==length(du) "Length of `u` is not equal to length of `du`."
+    u, t = munge_data(u, t)
+    p = CubicHermiteParameterCache(du, u, t)
+    return CubicHermiteSpline(du, u, t, p, extrapolate)
+end
+
+"""
+    QuinticHermiteSpline(ddu, du, u, t; extrapolate = false)
+
+It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polynomial such that the value and the first and second derivative are equal to given values in the data points.
+
+## Arguments
+
+  - `ddu`: the second derivative at the data points.
+  - `du`: the derivative at the data points.
+  - `u`: data points.
+  - `t`: time points.
+
+## Keyword Arguments
+
+  - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+"""
+struct QuinticHermiteSpline{uType, tType, duType, dduType, pType, T} <:
+       AbstractInterpolation{T}
+    ddu::uType
+    du::uType
+    u::uType
+    t::tType
+    p::QuinticHermiteParameterCache{pType}
+    extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
+    function QuinticHermiteSpline(ddu, du, u, t, p, extrapolate)
+        new{typeof(u), typeof(t), typeof(du), typeof(ddu), typeof(p.c₁), eltype(u)}(
+            ddu, du, u, t, p, extrapolate, Ref(1))
+    end
+end
+
+function QuinticHermiteSpline(ddu, du, u, t; extrapolate = false)
+    @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)
+    p = QuinticHermiteParameterCache(ddu, du, u, t)
+    return QuinticHermiteSpline(ddu, du, u, t, p, extrapolate)
+end
diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl
index 85344c93..0a543d8f 100644
--- a/src/interpolation_methods.jl
+++ b/src/interpolation_methods.jl
@@ -201,3 +201,25 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i
     end
     ucum, idx
 end
+
+# Cubic Hermite Spline
+function _interpolate(
+        A::CubicHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess)
+    Δt₀ = t - A.t[idx]
+    Δt₁ = t - A.t[idx + 1]
+    out = A.u[idx] + Δt₀ * A.du[idx]
+    out += Δt₀^2 * (A.p.c₁[idx] + Δt₁ * A.p.c₂[idx])
+    out, idx
+end
+
+# Quintic Hermite Spline
+function _interpolate(
+        A::QuinticHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess)
+    Δt₀ = t - A.t[idx]
+    Δt₁ = t - A.t[idx + 1]
+    out = A.u[idx] + Δt₀ * (A.du[idx] + A.ddu[idx] * Δt₀ / 2)
+    out += Δt₀^3 * (A.p.c₁[idx] + Δt₁ * (A.p.c₂[idx] + A.p.c₃[idx] * Δt₁))
+    out, idx
+end
diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl
new file mode 100644
index 00000000..66aa218f
--- /dev/null
+++ b/src/parameter_caches.jl
@@ -0,0 +1,49 @@
+struct CubicHermiteParameterCache{pType}
+    c₁::pType
+    c₂::pType
+end
+
+function CubicHermiteParameterCache(du, u, t)
+    parameters = CubicHermiteSplineParameters.(
+        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)
+    Δt = t[idx + 1] - t[idx]
+    u₀ = u[idx]
+    u₁ = u[idx + 1]
+    du₀ = du[idx]
+    du₁ = du[idx + 1]
+    c₁ = (u₁ - u₀ - du₀ * Δt) / Δt^2
+    c₂ = (du₁ - du₀ - 2c₁ * Δt) / Δt^2
+    return c₁, c₂
+end
+
+struct QuinticHermiteParameterCache{pType}
+    c₁::pType
+    c₂::pType
+    c₃::pType
+end
+
+function QuinticHermiteParameterCache(ddu, du, u, t)
+    parameters = QuinticHermiteSplineParameters.(
+        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)
+    Δt = t[idx + 1] - t[idx]
+    u₀ = u[idx]
+    u₁ = u[idx + 1]
+    du₀ = du[idx]
+    du₁ = du[idx + 1]
+    ddu₀ = ddu[idx]
+    ddu₁ = ddu[idx + 1]
+    c₁ = (u₁ - u₀ - du₀ * Δt - ddu₀ * Δt^2 / 2) / Δt^3
+    c₂ = (3u₀ - 3u₁ + 2(du₀ + du₁ / 2)Δt + ddu₀ * Δt^2 / 2) / Δt^4
+    c₃ = (6u₁ - 6u₀ - 3(du₀ + du₁)Δt + (ddu₁ - ddu₀)Δt^2 / 2) / Δt^5
+    return c₁, c₂, c₃
+end
diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl
index d3abbdd2..225ad9f3 100644
--- a/test/derivative_tests.jl
+++ b/test/derivative_tests.jl
@@ -8,8 +8,9 @@ using RegularizationTools
 using Optim
 using ForwardDiff
 
-function test_derivatives(method, u, t; args = [], kwargs = [], name::String)
-    func = method(u, t, args...; kwargs..., extrapolate = true)
+function test_derivatives(method; args = [], kwargs = [], name::String)
+    func = method(args...; kwargs..., extrapolate = true)
+    (; t) = func
     trange = collect(range(minimum(t) - 5.0, maximum(t) + 5.0, step = 0.1))
     trange_exclude = filter(x -> !in(x, t), trange)
     @testset "$name" begin
@@ -25,7 +26,8 @@ function test_derivatives(method, u, t; args = [], kwargs = [], name::String)
 
         # Interpolation time points
         for _t in t[2:(end - 1)]
-            if func isa BSplineInterpolation || func isa BSplineApprox
+            if func isa BSplineInterpolation || func isa BSplineApprox ||
+               func isa CubicHermiteSpline
                 fdiff = forward_fdm(5, 1; geom = true)(func, _t)
                 fdiff2 = forward_fdm(5, 1; geom = true)(t -> derivative(func, t), _t)
             else
@@ -65,7 +67,7 @@ function test_derivatives(method, u, t; args = [], kwargs = [], name::String)
     end
     @test_throws DataInterpolations.DerivativeNotFoundError derivative(
         func, t[1], 3)
-    func = method(u, t, args...)
+    func = method(args...)
     @test_throws DataInterpolations.ExtrapolationError derivative(func, t[1] - 1.0)
     @test_throws DataInterpolations.ExtrapolationError derivative(func, t[end] + 1.0)
     @test_throws DataInterpolations.DerivativeNotFoundError derivative(
@@ -75,52 +77,47 @@ end
 @testset "Linear Interpolation" begin
     u = vcat(collect(1:5), 2 * collect(6:10))
     t = 1.0collect(1:10)
-    test_derivatives(LinearInterpolation, u, t; name = "Linear Interpolation (Vector)")
+    test_derivatives(
+        LinearInterpolation; args = [u, t], name = "Linear Interpolation (Vector)")
     u = vcat(2.0collect(1:10)', 3.0collect(1:10)')
-    test_derivatives(LinearInterpolation, u, t; name = "Linear Interpolation (Matrix)")
+    test_derivatives(
+        LinearInterpolation; args = [u, t], name = "Linear Interpolation (Matrix)")
 end
 
 @testset "Quadratic Interpolation" begin
     u = [1.0, 4.0, 9.0, 16.0]
     t = [1.0, 2.0, 3.0, 4.0]
-    test_derivatives(QuadraticInterpolation,
-        u,
-        t;
+    test_derivatives(QuadraticInterpolation, args = [u, t],
         name = "Quadratic Interpolation (Vector)")
-    test_derivatives(QuadraticInterpolation,
-        u,
-        t;
-        args = [:Backward],
+    test_derivatives(QuadraticInterpolation;
+        args = [u, t, :Backward],
         name = "Quadratic Interpolation (Vector), backward")
     u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0]
-    test_derivatives(QuadraticInterpolation,
-        u,
-        t;
+    test_derivatives(QuadraticInterpolation;
+        args = [u, t],
         name = "Quadratic Interpolation (Matrix)")
 end
 
 @testset "Lagrange Interpolation" begin
     u = [1.0, 4.0, 9.0]
     t = [1.0, 2.0, 3.0]
-    test_derivatives(LagrangeInterpolation, u, t; name = "Lagrange Interpolation (Vector)")
+    test_derivatives(
+        LagrangeInterpolation; args = [u, t], name = "Lagrange Interpolation (Vector)")
     u = [1.0 4.0 9.0; 1.0 2.0 3.0]
-    test_derivatives(LagrangeInterpolation, u, t; name = "Lagrange Interpolation (Matrix)")
+    test_derivatives(
+        LagrangeInterpolation; args = [u, t], name = "Lagrange Interpolation (Matrix)")
     u = [[1.0, 4.0, 9.0], [3.0, 7.0, 4.0], [5.0, 4.0, 1.0]]
-    test_derivatives(LagrangeInterpolation,
-        u,
-        t;
+    test_derivatives(LagrangeInterpolation; args = [u, t],
         name = "Lagrange Interpolation (Vector of Vectors)")
     u = [[3.0 1.0 4.0; 1.0 5.0 9.0], [2.0 6.0 5.0; 3.0 5.0 8.0], [9.0 7.0 9.0; 3.0 2.0 3.0]]
-    test_derivatives(LagrangeInterpolation,
-        u,
-        t;
+    test_derivatives(LagrangeInterpolation; args = [u, t],
         name = "Lagrange Interpolation (Vector of Matrices)")
 end
 
 @testset "Akima 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)
-    test_derivatives(AkimaInterpolation, u, t; name = "Akima Interpolation")
+    test_derivatives(AkimaInterpolation; args = [u, t], name = "Akima Interpolation")
     @testset "Akima smooth derivative at end points" begin
         A = AkimaInterpolation(u, t)
         @test derivative(A, t[1]) ≈ derivative(A, nextfloat(t[1]))
@@ -131,56 +128,44 @@ end
 @testset "Quadratic Spline" begin
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
-    test_derivatives(QuadraticSpline, u, t; name = "Quadratic Interpolation (Vector)")
+    test_derivatives(
+        QuadraticSpline; args = [u, t], name = "Quadratic Interpolation (Vector)")
     u = [[1.0, 2.0, 9.0], [3.0, 7.0, 5.0], [5.0, 4.0, 1.0]]
-    test_derivatives(QuadraticSpline,
-        u,
-        t;
+    test_derivatives(QuadraticSpline; args = [u, t],
         name = "Quadratic Interpolation (Vector of Vectors)")
     u = [[1.0 4.0 9.0; 5.0 9.0 2.0], [3.0 7.0 4.0; 6.0 5.0 3.0], [5.0 4.0 1.0; 2.0 3.0 8.0]]
-    test_derivatives(QuadraticSpline,
-        u,
-        t;
+    test_derivatives(QuadraticSpline; args = [u, t],
         name = "Quadratic Interpolation (Vector of Matrices)")
 end
 
 @testset "Cubic Spline" begin
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
-    test_derivatives(CubicSpline, u, t; name = "Cubic Spline Interpolation (Vector)")
+    test_derivatives(
+        CubicSpline; args = [u, t], name = "Cubic Spline Interpolation (Vector)")
     u = [[1.0, 2.0, 9.0], [3.0, 7.0, 5.0], [5.0, 4.0, 1.0]]
-    test_derivatives(CubicSpline,
-        u,
-        t;
+    test_derivatives(CubicSpline; args = [u, t],
         name = "Cubic Spline Interpolation (Vector of Vectors)")
     u = [[1.0 4.0 9.0; 5.0 9.0 2.0], [3.0 7.0 4.0; 6.0 5.0 3.0], [5.0 4.0 1.0; 2.0 3.0 8.0]]
-    test_derivatives(CubicSpline,
-        u,
-        t;
+    test_derivatives(CubicSpline; args = [u, t],
         name = "Cubic Spline Interpolation (Vector of Matrices)")
 end
 
 @testset "BSplines" begin
     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]
-    test_derivatives(BSplineInterpolation,
-        u,
-        t;
-        args = [2,
+    test_derivatives(BSplineInterpolation;
+        args = [u, t, 2,
             :Uniform,
             :Uniform],
         name = "BSpline Interpolation (Uniform, Uniform)")
-    test_derivatives(BSplineInterpolation,
-        u,
-        t;
-        args = [2,
+    test_derivatives(BSplineInterpolation;
+        args = [u, t, 2,
             :ArcLen,
             :Average],
         name = "BSpline Interpolation (Arclen, Average)")
-    test_derivatives(BSplineApprox,
-        u,
-        t;
-        args = [
+    test_derivatives(BSplineApprox;
+        args = [u, t,
             3,
             4,
             :Uniform,
@@ -188,6 +173,32 @@ end
         name = "BSpline Approx (Uniform, Uniform)")
 end
 
+@testset "Cubic Hermite Spline" begin
+    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]
+    test_derivatives(CubicHermiteSpline; args = [du, u, t],
+        name = "Cubic Hermite Spline")
+    A = CubicHermiteSpline(du, u, t; extrapolate = true)
+    @test derivative.(Ref(A), t) ≈ du
+    @test derivative(A, 100.0)≈0.0105409 rtol=1e-5
+    @test derivative(A, 300.0)≈-0.0806717 rtol=1e-5
+end
+
+@testset "Quintic Hermite Spline" begin
+    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]
+    t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3]
+    test_derivatives(QuinticHermiteSpline; args = [ddu, du, u, t],
+        name = "Quintic Hermite Spline")
+    A = QuinticHermiteSpline(ddu, du, u, t; extrapolate = true)
+    @test derivative.(Ref(A), t) ≈ du
+    @test derivative.(Ref(A), t, 2) ≈ ddu
+    @test derivative(A, 100.0)≈0.0103916 rtol=1e-5
+    @test derivative(A, 300.0)≈0.0331361 rtol=1e-5
+end
+
 @testset "RegularizationSmooth" begin
     npts = 50
     xmin = 0.0
@@ -207,9 +218,7 @@ end
     tₒ = t[idx]
     uₒ = u[idx]
     A = RegularizationSmooth(uₒ, tₒ; alg = :fixed)
-    test_derivatives(RegularizationSmooth,
-        uₒ,
-        tₒ;
+    test_derivatives(RegularizationSmooth; args = [uₒ, tₒ],
         kwargs = [:alg => :fixed],
         name = "RegularizationSmooth")
 end
@@ -220,7 +229,7 @@ end
     t = range(-10, stop = 10, length = 40)
     u = model(t, [1.0, 2.0]) + 0.01 * randn(rng, length(t))
     p0 = [0.5, 0.5]
-    test_derivatives(Curvefit, u, t; args = [model, p0, LBFGS()], name = "Curvefit")
+    test_derivatives(Curvefit; args = [u, t, model, p0, LBFGS()], name = "Curvefit")
 end
 
 @testset "Symbolic derivatives" begin
diff --git a/test/integral_tests.jl b/test/integral_tests.jl
index 86a5c17f..a5641dff 100644
--- a/test/integral_tests.jl
+++ b/test/integral_tests.jl
@@ -5,8 +5,9 @@ using Optim, ForwardDiff
 using RegularizationTools
 using StableRNGs
 
-function test_integral(method, u, t; args = [], kwargs = [], name::String)
-    func = method(u, t, args...; kwargs..., extrapolate = true)
+function test_integral(method; args = [], kwargs = [], name::String)
+    func = method(args...; kwargs..., extrapolate = true)
+    (; t) = func
     t1 = minimum(t)
     t2 = maximum(t)
     @testset "$name" begin
@@ -46,7 +47,7 @@ function test_integral(method, u, t; args = [], kwargs = [], name::String)
         aint = integral(func, (t1 + t2) / 2, t2 + 5.0)
         @test isapprox(qint, aint, atol = 1e-6, rtol = 1e-8)
     end
-    func = method(u, t, args...; kwargs...)
+    func = method(args...; kwargs...)
     @test_throws DataInterpolations.ExtrapolationError integral(func, t[1] - 1.0)
     @test_throws DataInterpolations.ExtrapolationError integral(func, t[end] + 1.0)
     @test_throws DataInterpolations.ExtrapolationError integral(func, t[1] - 1.0, t[2])
@@ -56,27 +57,27 @@ end
 @testset "LinearInterpolation" begin
     u = 2.0collect(1:10)
     t = 1.0collect(1:10)
-    test_integral(LinearInterpolation, u, t; name = "Linear Interpolation (Vector)")
+    test_integral(
+        LinearInterpolation; args = [u, t], name = "Linear Interpolation (Vector)")
     u = round.(rand(100), digits = 5)
     t = 1.0collect(1:100)
-    test_integral(LinearInterpolation, u, t;
+    test_integral(LinearInterpolation; args = [u, t],
         name = "Linear Interpolation (Vector) with random points")
 end
 
 @testset "QuadraticInterpolation" begin
     u = [1.0, 4.0, 9.0, 16.0]
     t = [1.0, 2.0, 3.0, 4.0]
-    test_integral(QuadraticInterpolation, u, t; name = "Quadratic Interpolation (Vector)")
+    test_integral(
+        QuadraticInterpolation; args = [u, t], name = "Quadratic Interpolation (Vector)")
     u = [3.0, 0.0, 3.0, 0.0]
     t = [1.0, 2.0, 3.0, 4.0]
-    test_integral(QuadraticInterpolation,
-        u,
-        t;
-        args = [:Backward],
+    test_integral(QuadraticInterpolation;
+        args = [u, t, :Backward],
         name = "Quadratic Interpolation (Vector)")
     u = round.(rand(100), digits = 5)
     t = 1.0collect(1:10)
-    test_integral(QuadraticInterpolation, u, t;
+    test_integral(QuadraticInterpolation; args = [u, t],
         name = "Quadratic Interpolation (Vector) with random points")
 end
 
@@ -91,29 +92,64 @@ end
 @testset "QuadraticSpline" begin
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
-    test_integral(QuadraticSpline, u, t; name = "Quadratic Spline (Vector)")
+    test_integral(QuadraticSpline; args = [u, t], name = "Quadratic Spline (Vector)")
     u = round.(rand(100), digits = 5)
     t = 1.0collect(1:100)
     test_integral(
-        QuadraticSpline, u, t; name = "Quadratic Spline (Vector) with random points")
+        QuadraticSpline; args = [u, t], name = "Quadratic Spline (Vector) with random points")
 end
 
 @testset "CubicSpline" begin
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
-    test_integral(CubicSpline, u, t; name = "Cubic Spline (Vector)")
+    test_integral(CubicSpline; args = [u, t], name = "Cubic Spline (Vector)")
     u = round.(rand(100), digits = 5)
     t = 1.0collect(1:100)
-    test_integral(CubicSpline, u, t; name = "Cubic Spline (Vector) with random points")
+    test_integral(
+        CubicSpline; args = [u, t], name = "Cubic Spline (Vector) with random points")
 end
 
 @testset "AkimaInterpolation" 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)
-    test_integral(AkimaInterpolation, u, t; name = "Akima Interpolation (Vector)")
+    test_integral(AkimaInterpolation; args = [u, t], name = "Akima Interpolation (Vector)")
+    u = round.(rand(100), digits = 5)
+    t = 1.0collect(1:100)
+    test_integral(
+        AkimaInterpolation; args = [u, t], name = "Akima Interpolation (Vector) with random points")
+end
+
+@testset "CubicHermiteSpline" begin
+    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]
+    test_integral(CubicHermiteSpline; args = [du, u, t],
+        name = "Cubic Hermite Spline (Vector)")
+
+    u = round.(rand(100), digits = 5)
+    t = 1.0collect(1:100)
+    du = diff(u) ./ diff(t)
+    push!(du, 0)
+    test_integral(CubicHermiteSpline; args = [du, u, t],
+        name = "Cubic Hermite Spline (Vector) with random points")
+end
+
+@testset "QuinticHermiteSpline" begin
+    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]
+    t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3]
+    test_integral(QuinticHermiteSpline; args = [ddu, du, u, t],
+        name = "Quintic Hermite Spline (Vector)")
+
     u = round.(rand(100), digits = 5)
     t = 1.0collect(1:100)
-    test_integral(AkimaInterpolation, u, t; name = "Akima Interpolation (Vector)")
+    du = diff(u) ./ diff(t)
+    push!(du, 0)
+    ddu = diff(du) ./ diff(t)
+    push!(ddu, 0)
+    test_integral(QuinticHermiteSpline; args = [ddu, du, u, t],
+        name = "Quintic Hermite Spline (Vector) with random points")
 end
 
 @testset "RegularizationSmooth" begin
@@ -134,9 +170,8 @@ end
     idx = sortperm(t)
     tₒ = t[idx]
     uₒ = u[idx]
-    test_integral(RegularizationSmooth,
-        uₒ,
-        tₒ;
+    test_integral(RegularizationSmooth;
+        args = [uₒ, tₒ],
         kwargs = [:alg => :fixed],
         name = "RegularizationSmooth")
 end
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index a7c4eb2e..f9ef87f8 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -622,6 +622,33 @@ end
     end
 end
 
+@testset "Cubic Hermite Spline" begin
+    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]
+    A = CubicHermiteSpline(du, u, t; extrapolate = true)
+    @test A.(t) ≈ u
+    @test A(100.0)≈10.106770 rtol=1e-5
+    @test A(300.0)≈9.901542 rtol=1e-5
+    test_cached_index(A)
+    push!(u, 1.0)
+    @test_throws AssertionError CubicHermiteSpline(du, u, t)
+end
+
+@testset "Quintic Hermite Spline" begin
+    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]
+    t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3]
+    A = QuinticHermiteSpline(ddu, du, u, t; extrapolate = true)
+    @test A.(t) ≈ u
+    @test A(100.0)≈10.107996 rtol=1e-5
+    @test A(300.0)≈11.364162 rtol=1e-5
+    test_cached_index(A)
+    push!(u, 1.0)
+    @test_throws AssertionError QuinticHermiteSpline(ddu, du, u, t)
+end
+
 @testset "Curvefit" begin
     # Curvefit Interpolation
     rng = StableRNG(12345)
diff --git a/test/parameter_tests.jl b/test/parameter_tests.jl
new file mode 100644
index 00000000..3bedf3bc
--- /dev/null
+++ b/test/parameter_tests.jl
@@ -0,0 +1,21 @@
+using DataInterpolations
+
+@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]
+    t = collect(1:5)
+    A = CubicHermiteSpline(du, u, t)
+    @test A.p.c₁ ≈ [-1.0, -5.0, -5.0, -8.0]
+    @test A.p.c₂ ≈ [0.0, 13.0, 12.0, 9.0]
+end
+
+@testset "Quintic Hermite Spline" begin
+    ddu = [0.0, 3.0, 6.0, 4.0, 5.0]
+    du = [5.0, 3.0, 6.0, 8.0, 1.0]
+    u = [1.0, 5.0, 3.0, 4.0, 4.0]
+    t = collect(1:5)
+    A = QuinticHermiteSpline(ddu, du, u, t)
+    @test A.p.c₁ ≈ [-1.0, -6.5, -8.0, -10.0]
+    @test A.p.c₂ ≈ [1.0, 19.5, 20.0, 19.0]
+    @test A.p.c₃ ≈ [1.5, -37.5, -37.0, -26.5]
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index a7bac9e3..01f3a07e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,6 +2,7 @@ using SafeTestsets
 
 @safetestset "Quality Assurance" include("qa.jl")
 @safetestset "Interface" include("interface.jl")
+@safetestset "Parameter Tests" include("parameter_tests.jl")
 @safetestset "Interpolation Tests" include("interpolation_tests.jl")
 @safetestset "Derivative Tests" include("derivative_tests.jl")
 @safetestset "Integral Tests" include("integral_tests.jl")