diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl
index 5f035e62..393f1160 100644
--- a/src/DataInterpolations.jl
+++ b/src/DataInterpolations.jl
@@ -2,7 +2,7 @@ module DataInterpolations
 
 ### Interface Functionality
 
-abstract type AbstractInterpolation{T} end
+abstract type AbstractInterpolation{T, N} end
 
 using LinearAlgebra, RecipesBase
 using PrettyTables
@@ -92,8 +92,8 @@ export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation,
 
 # added for RegularizationSmooth, JJS 11/27/21
 ### Regularization data smoothing and interpolation
-struct RegularizationSmooth{uType, tType, T, T2, ITP <: AbstractInterpolation{T}} <:
-       AbstractInterpolation{T}
+struct RegularizationSmooth{uType, tType, T, T2, N, ITP <: AbstractInterpolation{T, N}} <:
+       AbstractInterpolation{T, N}
     u::uType
     û::uType
     t::tType
@@ -116,7 +116,8 @@ struct RegularizationSmooth{uType, tType, T, T2, ITP <: AbstractInterpolation{T}
             alg,
             Aitp,
             extrapolate)
-        new{typeof(u), typeof(t), eltype(u), typeof(λ), typeof(Aitp)}(
+        N = get_output_dim(u)
+        new{typeof(u), typeof(t), eltype(u), typeof(λ), N, typeof(Aitp)}(
             u,
             û,
             t,
@@ -143,8 +144,9 @@ struct CurvefitCache{
     lbType,
     algType,
     pminType,
-    T
-} <: AbstractInterpolation{T}
+    T,
+    N
+} <: AbstractInterpolation{T, N}
     u::uType
     t::tType
     m::mType        # model type
@@ -155,9 +157,10 @@ struct CurvefitCache{
     pmin::pminType  # optimized params
     extrapolate::Bool
     function CurvefitCache(u, t, m, p0, ub, lb, alg, pmin, extrapolate)
+        N = get_output_dim(u)
         new{typeof(u), typeof(t), typeof(m),
             typeof(p0), typeof(ub), typeof(lb),
-            typeof(alg), typeof(pmin), eltype(u)}(u,
+            typeof(alg), typeof(pmin), eltype(u), N}(u,
             t,
             m,
             p0,
diff --git a/src/derivatives.jl b/src/derivatives.jl
index bdb4d770..86640360 100644
--- a/src/derivatives.jl
+++ b/src/derivatives.jl
@@ -151,14 +151,14 @@ 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 = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.N
-    spline_coefficients!(N, A.d - 1, A.k, t_)
+    sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
+    spline_coefficients!(sc, 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])
     else
         for i in 1:(n - 1)
-            ducum += N[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
+            ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
         end
     end
     ducum * A.d * scale
@@ -172,14 +172,14 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig
     idx = get_idx(A, 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 = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.N
-    spline_coefficients!(N, A.d - 1, A.k, t_)
+    sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
+    spline_coefficients!(sc, 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])
     else
         for i in 1:(A.h - 1)
-            ducum += N[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
+            ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
         end
     end
     ducum * A.d * scale
diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl
index 4cddf7d3..33621f1a 100644
--- a/src/integral_inverses.jl
+++ b/src/integral_inverses.jl
@@ -1,4 +1,4 @@
-abstract type AbstractIntegralInverseInterpolation{T} <: AbstractInterpolation{T} end
+abstract type AbstractIntegralInverseInterpolation{T, N} <: AbstractInterpolation{T, N} end
 
 """
     invert_integral(A::AbstractInterpolation)::AbstractIntegralInverseInterpolation
@@ -33,15 +33,16 @@ Can be easily constructed with `invert_integral(A::LinearInterpolation{<:Abstrac
   - `t` : Given by `A.I` (the cumulative integral of `A`)
   - `A` : The `LinearInterpolation` object
 """
-struct LinearInterpolationIntInv{uType, tType, itpType, T} <:
-       AbstractIntegralInverseInterpolation{T}
+struct LinearInterpolationIntInv{uType, tType, itpType, T, N} <:
+       AbstractIntegralInverseInterpolation{T, N}
     u::uType
     t::tType
     extrapolate::Bool
     iguesser::Guesser{tType}
     itp::itpType
     function LinearInterpolationIntInv(u, t, A)
-        new{typeof(u), typeof(t), typeof(A), eltype(u)}(
+        N = get_output_dim(u)
+        new{typeof(u), typeof(t), typeof(A), eltype(u), N}(
             u, t, A.extrapolate, Guesser(t), A)
     end
 end
@@ -79,15 +80,16 @@ Can be easily constructed with `invert_integral(A::ConstantInterpolation{<:Abstr
   - `t` : Given by `A.I` (the cumulative integral of `A`)
   - `A` : The `ConstantInterpolation` object
 """
-struct ConstantInterpolationIntInv{uType, tType, itpType, T} <:
-       AbstractIntegralInverseInterpolation{T}
+struct ConstantInterpolationIntInv{uType, tType, itpType, T, N} <:
+       AbstractIntegralInverseInterpolation{T, N}
     u::uType
     t::tType
     extrapolate::Bool
     iguesser::Guesser{tType}
     itp::itpType
     function ConstantInterpolationIntInv(u, t, A)
-        new{typeof(u), typeof(t), typeof(A), eltype(u)}(
+        N = get_output_dim(u)
+        new{typeof(u), typeof(t), typeof(A), eltype(u), N}(
             u, t, A.extrapolate, Guesser(t), A
         )
     end
diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index d7dce336..6d124dcf 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -19,7 +19,7 @@ Extrapolation extends the last linear polynomial on each side.
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T}
+struct LinearInterpolation{uType, tType, IType, pType, T, N} <: AbstractInterpolation{T, N}
     u::uType
     t::tType
     I::IType
@@ -30,7 +30,8 @@ struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolati
     linear_lookup::Bool
     function LinearInterpolation(u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
-        new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}(
+        N = get_output_dim(u)
+        new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u), N}(
             u, t, I, p, extrapolate, Guesser(t), cache_parameters, linear_lookup)
     end
 end
@@ -66,7 +67,8 @@ Extrapolation extends the last quadratic polynomial on each side.
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T}
+struct QuadraticInterpolation{uType, tType, IType, pType, T, N} <:
+       AbstractInterpolation{T, N}
     u::uType
     t::tType
     I::IType
@@ -81,7 +83,8 @@ struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpol
         mode ∈ (:Forward, :Backward) ||
             error("mode should be :Forward or :Backward for QuadraticInterpolation")
         linear_lookup = seems_linear(assume_linear_t, t)
-        new{typeof(u), typeof(t), typeof(I), typeof(p.l₀), eltype(u)}(
+        N = get_output_dim(u)
+        new{typeof(u), typeof(t), typeof(I), typeof(p.l₀), eltype(u), N}(
             u, t, I, p, mode, extrapolate, Guesser(t), cache_parameters, linear_lookup)
     end
 end
@@ -116,8 +119,8 @@ It is the method of interpolation using Lagrange polynomials of (k-1)th order pa
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
 """
-struct LagrangeInterpolation{uType, tType, T, bcacheType} <:
-       AbstractInterpolation{T}
+struct LagrangeInterpolation{uType, tType, T, bcacheType, N} <:
+       AbstractInterpolation{T, N}
     u::uType
     t::tType
     n::Int
@@ -129,7 +132,8 @@ struct LagrangeInterpolation{uType, tType, T, bcacheType} <:
         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,
+        N = get_output_dim(u)
+        new{typeof(u), typeof(t), eltype(u), typeof(bcache), N}(u,
             t,
             n,
             bcache,
@@ -169,8 +173,8 @@ Extrapolation extends the last cubic polynomial on each side.
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
-       AbstractInterpolation{T}
+struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T, N} <:
+       AbstractInterpolation{T, N}
     u::uType
     t::tType
     I::IType
@@ -184,8 +188,9 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
     function AkimaInterpolation(
             u, t, I, b, c, d, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
+        N = get_output_dim(u)
         new{typeof(u), typeof(t), typeof(I), typeof(b), typeof(c),
-            typeof(d), eltype(u)}(u,
+            typeof(d), eltype(u), N}(u,
             t,
             I,
             b,
@@ -251,7 +256,7 @@ Extrapolation extends the last constant polynomial at the end points on each sid
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T}
+struct ConstantInterpolation{uType, tType, IType, T, N} <: AbstractInterpolation{T, N}
     u::uType
     t::tType
     I::IType
@@ -264,7 +269,8 @@ struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T}
     function ConstantInterpolation(
             u, t, I, dir, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
-        new{typeof(u), typeof(t), typeof(I), eltype(u)}(
+        N = get_output_dim(u)
+        new{typeof(u), typeof(t), typeof(I), eltype(u), N}(
             u, t, I, nothing, dir, extrapolate, Guesser(t), cache_parameters, linear_lookup)
     end
 end
@@ -299,8 +305,8 @@ Extrapolation extends the last quadratic polynomial on each side.
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <:
-       AbstractInterpolation{T}
+struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T, N} <:
+       AbstractInterpolation{T, N}
     u::uType
     t::tType
     I::IType
@@ -315,8 +321,9 @@ struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <:
     function QuadraticSpline(
             u, t, I, p, tA, d, z, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
+        N = get_output_dim(u)
         new{typeof(u), typeof(t), typeof(I), typeof(p.σ), typeof(tA),
-            typeof(d), typeof(z), eltype(u)}(u,
+            typeof(d), typeof(z), eltype(u), N}(u,
             t,
             I,
             p,
@@ -402,7 +409,8 @@ Second derivative on both ends are zero, which are also called "natural" boundar
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInterpolation{T}
+struct CubicSpline{uType, tType, IType, pType, hType, zType, T, N} <:
+       AbstractInterpolation{T, N}
     u::uType
     t::tType
     I::IType
@@ -415,7 +423,9 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInter
     linear_lookup::Bool
     function CubicSpline(u, t, I, p, h, z, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
-        new{typeof(u), typeof(t), typeof(I), typeof(p.c₁), typeof(h), typeof(z), eltype(u)}(
+        N = get_output_dim(u)
+        new{typeof(u), typeof(t), typeof(I), typeof(p.c₁),
+            typeof(h), typeof(z), eltype(u), N}(
             u,
             t,
             I,
@@ -508,15 +518,15 @@ Extrapolation is a constant polynomial of the end points on each side.
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
-       AbstractInterpolation{T}
+struct BSplineInterpolation{uType, tType, pType, kType, cType, scType, T, N} <:
+       AbstractInterpolation{T, N}
     u::uType
     t::tType
     d::Int    # degree
     p::pType  # params vector
     k::kType  # knot vector
     c::cType  # control points
-    N::NType  # Spline coefficients (preallocated memory)
+    sc::scType  # Spline coefficients (preallocated memory)
     pVecType::Symbol
     knotVecType::Symbol
     extrapolate::Bool
@@ -528,19 +538,21 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
             p,
             k,
             c,
-            N,
+            sc,
             pVecType,
             knotVecType,
             extrapolate,
             assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
-        new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(N), eltype(u)}(u,
+        N = get_output_dim(u)
+        new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(sc), eltype(u), N}(
+            u,
             t,
             d,
             p,
             k,
             c,
-            N,
+            sc,
             pVecType,
             knotVecType,
             extrapolate,
@@ -611,12 +623,12 @@ function BSplineInterpolation(
         end
     end
     # control points
-    N = zeros(eltype(t), n, n)
-    spline_coefficients!(N, d, k, p)
-    c = vec(N \ u[:, :])
-    N = zeros(eltype(t), n)
+    sc = zeros(eltype(t), n, n)
+    spline_coefficients!(sc, d, k, p)
+    c = vec(sc \ u[:, :])
+    sc = zeros(eltype(t), n)
     BSplineInterpolation(
-        u, t, d, p, k, c, N, pVecType, knotVecType, extrapolate, assume_linear_t)
+        u, t, d, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t)
 end
 
 """
@@ -643,8 +655,8 @@ Extrapolation is a constant polynomial of the end points on each side.
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
-       AbstractInterpolation{T}
+struct BSplineApprox{uType, tType, pType, kType, cType, scType, T, N} <:
+       AbstractInterpolation{T, N}
     u::uType
     t::tType
     d::Int    # degree
@@ -652,7 +664,7 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
     p::pType  # params vector
     k::kType  # knot vector
     c::cType  # control points
-    N::NType  # Spline coefficients (preallocated memory)
+    sc::scType  # Spline coefficients (preallocated memory)
     pVecType::Symbol
     knotVecType::Symbol
     extrapolate::Bool
@@ -665,21 +677,23 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
             p,
             k,
             c,
-            N,
+            sc,
             pVecType,
             knotVecType,
             extrapolate,
             assume_linear_t
     )
         linear_lookup = seems_linear(assume_linear_t, t)
-        new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(N), eltype(u)}(u,
+        N = get_output_dim(u)
+        new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(sc), eltype(u), N}(
+            u,
             t,
             d,
             h,
             p,
             k,
             c,
-            N,
+            sc,
             pVecType,
             knotVecType,
             extrapolate,
@@ -755,28 +769,28 @@ function BSplineApprox(
     c[1] = u[1]
     c[end] = u[end]
     q = zeros(eltype(u), n)
-    N = zeros(eltype(t), n, h)
+    sc = zeros(eltype(t), n, h)
     for i in 1:n
-        spline_coefficients!(view(N, i, :), d, k, p[i])
+        spline_coefficients!(view(sc, 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]
+        q[k] = u[k] - sc[k, 1] * u[1] - sc[k, h] * u[end]
     end
     Q = Matrix{eltype(u)}(undef, h - 2, 1)
     for i in 2:(h - 1)
         s = 0.0
         for k in 2:(n - 1)
-            s += N[k, i] * q[k]
+            s += sc[k, i] * q[k]
         end
         Q[i - 1] = s
     end
-    N = N[2:(end - 1), 2:(h - 1)]
-    M = transpose(N) * N
+    sc = sc[2:(end - 1), 2:(h - 1)]
+    M = transpose(sc) * sc
     P = M \ Q
     c[2:(end - 1)] .= vec(P)
-    N = zeros(eltype(t), h)
+    sc = zeros(eltype(t), h)
     BSplineApprox(
-        u, t, d, h, p, k, c, N, pVecType, knotVecType, extrapolate, assume_linear_t)
+        u, t, d, h, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t)
 end
 
 """
@@ -799,7 +813,8 @@ It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomi
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInterpolation{T}
+struct CubicHermiteSpline{uType, tType, IType, duType, pType, T, N} <:
+       AbstractInterpolation{T, N}
     du::duType
     u::uType
     t::tType
@@ -812,7 +827,8 @@ struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInte
     function CubicHermiteSpline(
             du, u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
-        new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u)}(
+        N = get_output_dim(u)
+        new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u), N}(
             du, u, t, I, p, extrapolate, Guesser(t), cache_parameters, linear_lookup)
     end
 end
@@ -878,8 +894,8 @@ It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polyno
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <:
-       AbstractInterpolation{T}
+struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T, N} <:
+       AbstractInterpolation{T, N}
     ddu::dduType
     du::duType
     u::uType
@@ -893,8 +909,9 @@ struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <:
     function QuinticHermiteSpline(
             ddu, du, u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
+        N = get_output_dim(u)
         new{typeof(u), typeof(t), typeof(I), typeof(du),
-            typeof(ddu), typeof(p.c₁), eltype(u)}(
+            typeof(ddu), typeof(p.c₁), eltype(u), N}(
             ddu, du, u, t, I, p, extrapolate, Guesser(t), cache_parameters, linear_lookup)
     end
 end
diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl
index 03adf558..194c6f5b 100644
--- a/src/interpolation_methods.jl
+++ b/src/interpolation_methods.jl
@@ -176,11 +176,11 @@ function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}},
     idx = get_idx(A, 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 = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.N
-    nonzero_coefficient_idxs = spline_coefficients!(N, A.d, A.k, t)
+    sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
+    nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t)
     ucum = zero(eltype(A.u))
     for i in nonzero_coefficient_idxs
-        ucum += N[i] * A.c[i]
+        ucum += sc[i] * A.c[i]
     end
     ucum
 end
@@ -192,11 +192,11 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i
     # change t into param [0 1]
     idx = get_idx(A, 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 = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.N
-    nonzero_coefficient_idxs = spline_coefficients!(N, A.d, A.k, t)
+    sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
+    nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t)
     ucum = zero(eltype(A.u))
     for i in nonzero_coefficient_idxs
-        ucum += N[i] * A.c[i]
+        ucum += sc[i] * A.c[i]
     end
     ucum
 end
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index f8fc7147..57302257 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -59,6 +59,19 @@ function spline_coefficients!(N, d, k, u::AbstractVector)
     return nothing
 end
 
+# Get Output Dimension for parameterizing AbstractInterpolations
+function get_output_dim(u::AbstractVector{<:Number})
+    return (1,)
+end
+
+function get_output_dim(u::AbstractVector)
+    return (length(first(u)),)
+end
+
+function get_output_dim(u::AbstractArray)
+    return size(u)[1:(end - 1)]
+end
+
 # helper function for data manipulation
 function munge_data(u::AbstractVector{<:Real}, t::AbstractVector{<:Real})
     return u, t