diff --git a/Project.toml b/Project.toml
index 2e0267a66..4bdc82f01 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,20 +9,27 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 Aqua = "0.8"
 ChainRulesCore = "1.24"
 ControlSystemsBase = "1.4"
+DataFrames = "1.7"
 DataInterpolations = "6"
+ForwardDiff = "0.10"
 DiffEqBase = "6.152"
 IfElse = "0.1"
 LinearAlgebra = "1.10"
-ModelingToolkit = "9.46.1"
+Optimization = "4"
+ModelingToolkit = "9.47"
 OrdinaryDiffEq = "6.87"
 OrdinaryDiffEqDefault = "1.1"
+PreallocationTools = "0.4.23"
 SafeTestsets = "0.1"
+SciMLStructures = "1.4.2"
+SymbolicIndexingInterface = "0.3.28"
 Symbolics = "6.14"
 Test = "1"
 julia = "1.10"
@@ -30,12 +37,17 @@ julia = "1.10"
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 OrdinaryDiffEqDefault = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226"
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "LinearAlgebra", "OrdinaryDiffEq", "OrdinaryDiffEqDefault", "SafeTestsets", "Test", "ControlSystemsBase", "DataInterpolations"]
+test = ["Aqua", "LinearAlgebra", "OrdinaryDiffEqDefault", "OrdinaryDiffEq", "Optimization", "SafeTestsets", "Test", "ControlSystemsBase", "DataFrames", "DataInterpolations", "SciMLStructures", "SymbolicIndexingInterface", "ForwardDiff"]
diff --git a/docs/Project.toml b/docs/Project.toml
index dffb4fce9..2c3701147 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,5 +1,7 @@
 [deps]
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
@@ -10,6 +12,8 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
 [compat]
 ControlSystemsBase = "1.1"
+DataFrames = "1.7"
+DataInterpolations = "6.4"
 DifferentialEquations = "7.6"
 Documenter = "1"
 IfElse = "0.1"
diff --git a/docs/src/tutorials/dc_motor_pi.md b/docs/src/tutorials/dc_motor_pi.md
index 5b5f6982e..b4e84f38b 100644
--- a/docs/src/tutorials/dc_motor_pi.md
+++ b/docs/src/tutorials/dc_motor_pi.md
@@ -77,7 +77,7 @@ so that it can be represented as a system of `ODEs` (ordinary differential equat
 
 ```@example dc_motor_pi
 sys = structural_simplify(model)
-prob = ODEProblem(sys, unknowns(sys) .=> 0.0, (0, 6.0))
+prob = ODEProblem(sys, [sys.L1.i => 0.0], (0, 6.0))
 sol = solve(prob)
 
 p1 = plot(sol.t, sol[sys.inertia.w], ylabel = "Angular Vel. in rad/s",
diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md
index 414a22eed..f48dd96b1 100644
--- a/docs/src/tutorials/input_component.md
+++ b/docs/src/tutorials/input_component.md
@@ -1,59 +1,171 @@
-# Running Models with Discrete Data
+# Building Models with Discrete Data, Interpolations, and Lookup Tables
 
-There are 3 ways to include data as part of a model.
+There are 4 ways to include data as part of a model.
 
- 1. using `ModelingToolkitStandardLibrary.Blocks.TimeVaryingFunction`
- 2. using a custom component with external data
- 3. using `ModelingToolkitStandardLibrary.Blocks.SampledData`
+ 1. using `ModelingToolkitStandardLibrary.Blocks.Interpolation`
+ 2. using `ModelingToolkitStandardLibrary.Blocks.ParametrizedInterpolation`
+ 3. using a custom component with external data (not recommended)
+ 4. using `ModelingToolkitStandardLibrary.Blocks.SampledData` (legacy)
 
 This tutorial demonstrate each case and explain the pros and cons of each.
 
-## `TimeVaryingFunction` Component
+## `Interpolation` Block
 
-The `ModelingToolkitStandardLibrary.Blocks.TimeVaryingFunction` component is easy to use and is performant.  However the data is locked to the `ODESystem` and can only be changed by building a new `ODESystem`.  Therefore, running a batch of data would not be efficient.  Below is an example of how to use the `TimeVaryingFunction` with `DataInterpolations` to build the function from sampled discrete data.
+The `ModelingToolkitStandardLibrary.Blocks.Interpolation` component is easy to use and is performant.
+It is similar to using callable parameters, but it provides a block interface with `RealInput` and `RealOutput` connectors.
+The `Interpolation` is compatible with interpolation types from `DataInterpolation`.
 
-```julia
+```@docs
+ModelingToolkitStandardLibrary.Blocks.Interpolation
+```
+
+Here is an example on how to use it. Let's consider a mass-spring-damper system, where
+we have an external force as an input. We then generate some example data in a `DataFrame`
+that would represent a measurement of the input. In a more realistic case, this `DataFrame`
+would be read from a file.
+
+```@example interpolation_block
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using ModelingToolkitStandardLibrary.Blocks
 using DataInterpolations
 using OrdinaryDiffEq
+using DataFrames
+using Plots
 
-function System(f; name)
-    src = TimeVaryingFunction(f)
+function MassSpringDamper(; name)
+    @named input = RealInput()
+    @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0
+    @parameters m=10 k=1000 d=1
 
-    vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0
+    eqs = [
+           f ~ input.u
+           ddx * 10 ~ k * x + d * dx + f
+           D(x) ~ dx
+           D(dx) ~ ddx]
+
+    ODESystem(eqs, t; name, systems = [input])
+end
+
+function MassSpringDamperSystem(data, time; name)
+    @named src = Interpolation(LinearInterpolation, data, time)
+    @named clk = ContinuousClock()
+    @named model = MassSpringDamper()
+
+    eqs = [
+        connect(src.input, clk.output)
+        connect(src.output, model.input)
+    ]
+
+    ODESystem(eqs, t, [], []; name, systems = [src, clk, model])
+end
+
+function generate_data()
+    dt = 4e-4
+    time = 0:dt:0.1
+    data = sin.(2 * pi * time * 100)
+
+    return DataFrame(; time, data)
+end
+
+df = generate_data() # example data
+
+@named system = MassSpringDamperSystem(df.data, df.time)
+sys = structural_simplify(system)
+prob = ODEProblem(sys, [], (0, df.time[end]))
+sol = solve(prob)
+plot(sol)
+```
+
+## `ParametrizedInterpolation` Block
+
+The `ModelingToolkitStandardLibrary.Blocks.ParametrizedInterpolation` component is similar to `Interpolation`, but as the name suggests, it is parametrized by the data, allowing one to change the underlying data without rebuilding the model as the data is represented via vector parameters.
+The main advantage of this block over the [`Interpolation`](@ref) one is that one can use it for optimization problems. Currently, this supports forward mode AD via ForwardDiff, but due to the increased flexibility of the types in the component, this is not as fast as the `Interpolation` block,
+so it is recommended to use only when the added flexibility is required.
+
+```@docs
+ModelingToolkitStandardLibrary.Blocks.ParametrizedInterpolation
+```
+
+Here is an example on how to use it
+
+```@example parametrized_interpolation
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+using ModelingToolkitStandardLibrary.Blocks
+using DataInterpolations
+using OrdinaryDiffEq
+using DataFrames
+using Plots
+
+function MassSpringDamper(; name)
+    @named input = RealInput()
+    vars = @variables f(t) x(t)=0 dx(t) [guess=0] ddx(t)
     pars = @parameters m=10 k=1000 d=1
 
-    eqs = [f ~ src.output.u
+    eqs = [
+           f ~ input.u
            ddx * 10 ~ k * x + d * dx + f
            D(x) ~ dx
            D(dx) ~ ddx]
 
-    ODESystem(eqs, t, vars, pars; systems = [src], name)
+    ODESystem(eqs, t, vars, pars; name, systems = [input])
 end
 
-dt = 4e-4
-time = 0:dt:0.1
-data = sin.(2 * pi * time * 100) # example data
+function MassSpringDamperSystem(data, time; name)
+    @named src = ParametrizedInterpolation(LinearInterpolation, data, time)
+    @named clk = ContinuousClock()
+    @named model = MassSpringDamper()
+
+    eqs = [
+        connect(model.input, src.output)
+        connect(src.input, clk.output)
+    ]
 
-f = LinearInterpolation(data, time)
+    ODESystem(eqs, t; name, systems = [src, clk, model])
+end
+
+function generate_data()
+    dt = 4e-4
+    time = 0:dt:0.1
+    data = sin.(2 * pi * time * 100)
+
+    return DataFrame(; time, data)
+end
 
-@named system = System(f)
+df = generate_data() # example data
+
+@named system = MassSpringDamperSystem(df.data, df.time)
 sys = structural_simplify(system)
-prob = ODEProblem(sys, [], (0, time[end]))
-sol = solve(prob, ImplicitEuler())
+prob = ODEProblem(sys, [], (0, df.time[end]))
+sol = solve(prob)
+plot(sol)
+```
+
+If we want to run a new data set, this requires only remaking the problem and solving again
+```@example parametrized_interpolation
+prob2 = remake(prob, p = [sys.src.data => ones(length(df.data))])
+sol2 = solve(prob2)
+plot(sol2)
 ```
 
-If we want to run a new data set, this requires building a new `LinearInterpolation` and `ODESystem` followed by running `structural_simplify`, all of which takes time.  Therefore, to run several pieces of data it's better to re-use an `ODESystem`.  The next couple methods will demonstrate how to do this.
+!!! note
+    Note that when changing the data, the length of the new data must be the same as the length of the original data.
 
 ## Custom Component with External Data
 
 The below code shows how to include data using a `Ref` and registered `get_sampled_data` function.  This example uses a very basic function which requires non-adaptive solving and sampled data.  As can be seen, the data can easily be set and changed before solving.
 
-```julia
+```@example custom_component_external_data
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+using ModelingToolkitStandardLibrary.Blocks
+using OrdinaryDiffEq
+
 const rdata = Ref{Vector{Float64}}()
 
+dt = 4e-4
+time = 0:dt:0.1
 # Data Sets
 data1 = sin.(2 * pi * time * 100)
 data2 = cos.(2 * pi * time * 50)
@@ -106,8 +218,13 @@ Additional code could be added to resolve this issue, for example by using a `Re
 To resolve the issues presented above, the `ModelingToolkitStandardLibrary.Blocks.SampledData` component can be used which allows for a resusable `ODESystem` and self contained data which ensures a solution which remains valid for it's lifetime.  Now it's possible to also parallelize the call to `solve()`.
 
 ```julia
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+using ModelingToolkitStandardLibrary.Blocks
+using OrdinaryDiffEq
+
 function System(; name)
-    src = SampledData(Float64)
+    src = SampledData(Float64, name=:src)
 
     vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0
     pars = @parameters m=10 k=1000 d=1
@@ -121,16 +238,22 @@ function System(; name)
 end
 
 @named system = System()
-sys = structural_simplify(system)
+sys = structural_simplify(system, split=false)
 s = complete(system)
-prob = ODEProblem(sys, [], (0, time[end]); tofloat = false)
+
+dt = 4e-4
+time = 0:dt:0.1
+data1 = sin.(2 * pi * time * 100)
+data2 = cos.(2 * pi * time * 50)
+
+prob = ODEProblem(sys, [], (0, time[end]); split=false, tofloat = false, use_union=true)
 defs = ModelingToolkit.defaults(sys)
 
 function get_prob(data)
     defs[s.src.buffer] = Parameter(data, dt)
     # ensure p is a uniform type of Vector{Parameter{Float64}} (converting from Vector{Any})
     p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat = false))
-    remake(prob; p)
+    remake(prob; p, build_initializeprob=false)
 end
 
 prob1 = get_prob(data1)
diff --git a/src/Blocks/Blocks.jl b/src/Blocks/Blocks.jl
index 3eceab637..adb5db6e5 100644
--- a/src/Blocks/Blocks.jl
+++ b/src/Blocks/Blocks.jl
@@ -16,7 +16,8 @@ export Log, Log10
 include("math.jl")
 
 export Constant, TimeVaryingFunction, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine,
-       Square, Triangular, Parameter, SampledData
+       Square, Triangular, Parameter, SampledData,
+       Interpolation, ParametrizedInterpolation
 include("sources.jl")
 
 export Limiter, DeadZone, SlewRateLimiter
diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl
index 5930be1e9..31a2a291f 100644
--- a/src/Blocks/sources.jl
+++ b/src/Blocks/sources.jl
@@ -1,5 +1,6 @@
 using DiffEqBase
 import ChainRulesCore
+using PreallocationTools
 
 # Define and register smooth functions
 # These are "smooth" aka differentiable and avoid Gibbs effect
@@ -730,3 +731,144 @@ end
 function SampledData(; name, buffer, sample_time, circular_buffer)
     SampledData(SampledDataType.vector_based; name, buffer, sample_time, circular_buffer)
 end
+
+"""
+    Interpolation(interp_type, u, x, args...; name)
+
+Represent function interpolation symbolically as a block component.
+By default interpolation types from [`DataInterpolations.jl`](https://github.com/SciML/DataInterpolations.jl) are supported,
+but in general any callable type that builds the interpolation object via `itp = interpolation_type(u, x, args...)` and calls
+the interpolation with `itp(t)` should work. This does not need to represent an interpolation, it can be any type that satisfies
+the interface, such as lookup tables.
+# Arguments:
+  - `interp_type`: the type of the interpolation. For `DataInterpolations`,
+these would be any of [the available interpolations](https://github.com/SciML/DataInterpolations.jl?tab=readme-ov-file#available-interpolations),
+such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`.
+  - `u`: the data used for interpolation. For `DataInterpolations` this will be an `AbstractVector`
+  - `x`: the values that each data points correspond to, usually the times corresponding to each value in `u`.
+  - `args`: any other arguments needed to build the interpolation
+# Keyword arguments:
+  - `name`: the name of the component
+
+# Parameters:
+  - `interpolator`: the symbolic representation of the interpolation object, callable as `interpolator(t)`
+
+# Connectors:
+  - `input`: a [`RealInput`](@ref) connector corresponding to the input variable
+  - `output`: a [`RealOutput`](@ref) connector corresponding to the interpolated value
+"""
+function Interpolation(interp_type, u, x, args...; name)
+    itp = interp_type(u, x, args...)
+    Interpolation(itp; name)
+end
+
+function Interpolation(itp; name)
+    @parameters (interpolator::typeof(itp))(..) = itp
+    @named input = RealInput()
+    @named output = RealOutput()
+
+    eqs = [output.u ~ interpolator(input.u)]
+
+    ODESystem(
+        eqs, t, [], [interpolator]; name, systems = [input, output])
+end
+
+"""
+    CachedInterpolation
+
+This callable struct caches the calls to an interpolation object via PreallocationTools.
+"""
+struct CachedInterpolation{T, I, U, X, C}
+    interpolation_type::I
+    prev_u::U
+    prev_x::X
+    cache::C
+
+    function CachedInterpolation(interpolation_type, u, x, args)
+        # we need to copy the inputs to avoid aliasing
+        prev_u = DiffCache(copy(u))
+        # Interpolation points can be a range, but we want to be able
+        # to update the cache if needed (and setindex! is not defined on ranges)
+        # with a view from MTKParameters, so we collect to get a vector
+        prev_x = DiffCache(collect(copy(x)))
+        cache = GeneralLazyBufferCache() do (u, x)
+            interpolation_type(get_tmp(prev_u, u), get_tmp(prev_x, x), args...)
+        end
+        T = typeof(cache[(get_tmp(prev_u, u), get_tmp(prev_x, x))])
+        I = typeof(interpolation_type)
+        U = typeof(prev_u)
+        X = typeof(prev_x)
+        C = typeof(cache)
+
+        new{T, I, U, X, C}(interpolation_type, prev_u, prev_x, cache)
+    end
+end
+
+function (f::CachedInterpolation{T})(u, x, args) where {T}
+    (; prev_u, prev_x, cache, interpolation_type) = f
+
+    interp = @inbounds if (u, x) ≠ (get_tmp(prev_u, u), get_tmp(prev_x, x))
+        get_tmp(prev_u, u) .= u
+        get_tmp(prev_x, x) .= x
+        cache.bufs[(u, x)] = interpolation_type(
+            get_tmp(prev_u, u), get_tmp(prev_x, x), args...)
+    else
+        cache[(u, x)]
+    end
+
+    return interp
+end
+
+Base.nameof(::CachedInterpolation) = :CachedInterpolation
+
+@register_symbolic (f::CachedInterpolation)(u::AbstractArray, x::AbstractArray, args::Tuple)
+
+"""
+    ParametrizedInterpolation(interp_type, u, x, args...; name, t = ModelingToolkit.t_nounits)
+
+Represent function interpolation symbolically as a block component, with the interpolation data represented parametrically.
+By default interpolation types from [`DataInterpolations.jl`](https://github.com/SciML/DataInterpolations.jl) are supported,
+but in general any callable type that builds the interpolation object via `itp = interpolation_type(u, x, args...)` and calls
+the interpolation with `itp(t)` should work. This does not need to represent an interpolation, it can be any type that satisfies
+the interface, such as lookup tables.
+# Arguments:
+  - `interp_type`: the type of the interpolation. For `DataInterpolations`,
+these would be any of [the available interpolations](https://github.com/SciML/DataInterpolations.jl?tab=readme-ov-file#available-interpolations),
+such as `LinearInterpolation`, `ConstantInterpolation` or `CubicSpline`.
+  - `u`: the data used for interpolation. For `DataInterpolations` this will be an `AbstractVector`
+  - `x`: the values that each data points correspond to, usually the times corresponding to each value in `u`.
+  - `args`: any other arguments beeded to build the interpolation
+# Keyword arguments:
+  - `name`: the name of the component
+
+# Parameters:
+  - `data`: the symbolic representation of the data passed at construction time via `u`.
+  - `ts`: the symbolic representation of times corresponding to the data passed at construction time via `x`.
+
+# Connectors:
+  - `input`: a [`RealInput`](@ref) connector corresponding to the independent variable
+  - `output`: a [`RealOutput`](@ref) connector corresponding to the interpolated value
+"""
+function ParametrizedInterpolation(
+        interp_type::T, u::AbstractVector, x::AbstractVector, args...;
+        name) where {T}
+    build_interpolation = CachedInterpolation(interp_type, u, x, args)
+
+    @parameters data[1:length(x)] = u
+    @parameters ts[1:length(x)] = x
+    @parameters interpolation_type::T=interp_type [tunable = false]
+    @parameters (interpolator::interp_type)(..)::eltype(u)
+
+    @named input = RealInput()
+    @named output = RealOutput()
+
+    eqs = [output.u ~ interpolator(input.u)]
+
+    ODESystem(eqs, ModelingToolkit.t_nounits, [],
+        [data, ts, interpolation_type, interpolator];
+        parameter_dependencies = [
+            interpolator ~ build_interpolation(data, ts, args)
+        ],
+        systems = [input, output],
+        name)
+end
diff --git a/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl b/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl
index ce5e26907..3403ee022 100644
--- a/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl
+++ b/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl
@@ -1,5 +1,5 @@
 """
-Library to model iso-thermal compressible liquid fluid flow 
+Library to model iso-thermal compressible liquid fluid flow
 """
 module IsothermalCompressible
 
diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl
index c073385b1..8598278b6 100644
--- a/test/Blocks/sources.jl
+++ b/test/Blocks/sources.jl
@@ -5,6 +5,12 @@ using ModelingToolkitStandardLibrary.Blocks: smooth_sin, smooth_cos, smooth_damp
                                              smooth_square, smooth_step, smooth_ramp,
                                              smooth_triangular, triangular, square
 using OrdinaryDiffEq: ReturnCode.Success
+using DataInterpolations
+using DataFrames
+using SymbolicIndexingInterface
+using SciMLStructures: SciMLStructures, Tunable
+using Optimization
+using ForwardDiff
 
 @testset "Constant" begin
     @named src = Constant(k = 2)
@@ -412,8 +418,6 @@ end
 end
 
 @testset "SampledData" begin
-    using DataInterpolations
-
     dt = 4e-4
     t_end = 10.0
     time = 0:dt:t_end
@@ -478,3 +482,137 @@ end
         @test sol[ddy][end]≈2 atol=1e-3
     end
 end
+
+@testset "Interpolation" begin
+    @variables y(t) = 0
+    u = rand(15)
+    x = 0:14.0
+
+    @named i = Interpolation(LinearInterpolation, u, x)
+    eqs = [i.input.u ~ t, D(y) ~ i.output.u]
+
+    @named model = ODESystem(eqs, t, systems = [i])
+    sys = structural_simplify(model)
+
+    prob = ODEProblem{true, SciMLBase.FullSpecialize}(sys, [], (0.0, 4))
+    sol = solve(prob, Tsit5())
+
+    @test SciMLBase.successful_retcode(sol)
+end
+
+@testset "ParametrizedInterpolation" begin
+    @variables y(t) = 0
+    u = rand(15)
+    x = 0:14.0
+
+    @testset "LinearInterpolation" begin
+        @named i = ParametrizedInterpolation(LinearInterpolation, u, x)
+        eqs = [i.input.u ~ t, D(y) ~ i.output.u]
+
+        @named model = ODESystem(eqs, t, systems = [i])
+        sys = structural_simplify(model)
+
+        prob = ODEProblem{true, SciMLBase.FullSpecialize}(sys, [], (0.0, 4))
+        sol = solve(prob, Tsit5())
+
+        @test SciMLBase.successful_retcode(sol)
+
+        prob2 = remake(prob, p = [i.data => ones(15)])
+        sol2 = solve(prob2)
+
+        @test SciMLBase.successful_retcode(sol2)
+        @test all(only.(sol2.u) .≈ sol2.t) # the solution for y' = 1 is y(t) = t
+
+        set_data! = setp(prob2, i.data)
+        set_data!(prob2, zeros(15))
+        sol3 = solve(prob2)
+        @test SciMLBase.successful_retcode(sol3)
+        @test iszero(sol3)
+
+        function loss(x, p)
+            prob0, set_data! = p
+            ps = parameter_values(prob0)
+            arr, repack, alias = SciMLStructures.canonicalize(Tunable(), ps)
+            T = promote_type(eltype(x), eltype(arr))
+            promoted_ps = SciMLStructures.replace(Tunable(), ps, T.(arr))
+            prob = remake(prob0; p = promoted_ps)
+
+            set_data!(prob, x)
+            sol = solve(prob)
+            sum(abs2.(only.(sol.u) .- sol.t))
+        end
+
+        set_data! = setp(prob, i.data)
+        of = OptimizationFunction(loss, AutoForwardDiff())
+        op = OptimizationProblem(
+            of, u, (prob, set_data!), lb = zeros(15), ub = fill(2.0, 15))
+
+        # check that type changing works
+        @test length(ForwardDiff.gradient(x -> of(x, (prob, set_data!)), u)) == 15
+
+        r = solve(op, Optimization.LBFGS(), maxiters = 1000)
+        @test of(r.u, (prob, set_data!)) < of(u, (prob, set_data!))
+    end
+
+    @testset "BSplineInterpolation" begin
+        @named i = ParametrizedInterpolation(
+            BSplineInterpolation, u, x, 3, :Uniform, :Uniform)
+        eqs = [i.input.u ~ t, D(y) ~ i.output.u]
+
+        @named model = ODESystem(eqs, t, systems = [i])
+        sys = structural_simplify(model)
+
+        prob = ODEProblem(sys, [], (0.0, 4))
+        sol = solve(prob)
+
+        @test SciMLBase.successful_retcode(sol)
+    end
+
+    @testset "Initialization" begin
+        function MassSpringDamper(; name)
+            @named input = RealInput()
+            vars = @variables f(t) x(t)=0 dx(t) [guess = 0] ddx(t)
+            pars = @parameters m=10 k=1000 d=1
+
+            eqs = [f ~ input.u
+                   ddx * 10 ~ k * x + d * dx + f
+                   D(x) ~ dx
+                   D(dx) ~ ddx]
+
+            ODESystem(eqs, t, vars, pars; name, systems = [input])
+        end
+
+        function MassSpringDamperSystem(data, time; name)
+            @named src = ParametrizedInterpolation(LinearInterpolation, data, time)
+            @named clk = ContinuousClock()
+            @named model = MassSpringDamper()
+
+            eqs = [connect(model.input, src.output)
+                   connect(src.input, clk.output)]
+
+            ODESystem(eqs, t; name, systems = [src, clk, model])
+        end
+
+        function generate_data()
+            dt = 4e-4
+            time = 0:dt:0.1
+            data = sin.(2 * pi * time * 100)
+
+            return DataFrame(; time, data)
+        end
+
+        df = generate_data() # example data
+
+        @named system = MassSpringDamperSystem(df.data, df.time)
+        sys = structural_simplify(system)
+        prob = ODEProblem(sys, [], (0, df.time[end]))
+        sol = solve(prob)
+
+        @test SciMLBase.successful_retcode(sol)
+
+        prob2 = remake(prob, p = [sys.src.data => ones(length(df.data))])
+        sol2 = solve(prob2)
+
+        @test SciMLBase.successful_retcode(sol2)
+    end
+end