diff --git a/Project.toml b/Project.toml
index 25c88dc..a07b72b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["Huy D. Vo <vdhuy91@gmail.com> and contributors"]
 version = "0.1.3"
 
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
@@ -13,24 +14,14 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ReusePatterns = "a39b5e78-89b5-562b-97d8-70689129df0c"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
-Catalyst = "10, 11, 12"
-DataStructures = "0.18, 0.19"
-DifferentialEquations = "7"
-DocStringExtensions = "0.8"
-ForwardDiff = "0.10, 0.11"
-LinearAlgebra = "1"
-ModelingToolkit = "8"
-ReusePatterns = "0.2, 0.3"
-StaticArrays = "1"
-Sundials = "4"
-Symbolics = "4"
-julia = "^1.6"
+
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/examples/hog1p.jl b/examples/hog1p.jl
index e93b60e..580afa9 100644
--- a/examples/hog1p.jl
+++ b/examples/hog1p.jl
@@ -13,7 +13,7 @@ function Hog1p(t)
     signal = Ahog * (u / (1.0 + u / Mhog))^η
     return signal 
 end
-@parameters k01, k10, a, k12, k21, k23, k32, λ0, λ1, λ2, λ3, ktrans, γnuc, γcyt
+
 rn = @reaction_network begin 
     k01, G0 --> G1
     max(0, k10 - a*Hog1p(t)), G1 --> G0 
@@ -28,7 +28,7 @@ rn = @reaction_network begin
     γnuc, RNAnuc --> ∅
     ktrans, RNAnuc --> RNAcyt 
     γcyt, RNAcyt --> ∅
-end k01 k10 a k12 k21 k23 k32 λ0 λ1 λ2 λ3 γnuc ktrans γcyt
+end 
 
 param_values = Dict([
 k01=> 2.6e-3,
diff --git a/examples/telegraph_cme.jl b/examples/telegraph_cme.jl
index 80d3a0b..013b741 100644
--- a/examples/telegraph_cme.jl
+++ b/examples/telegraph_cme.jl
@@ -42,13 +42,12 @@ tspan = (0.0, 300.0)
 fspsol1 = solve(model, p0, tspan, fspalgorithm)
 
 # Bursting model definition using Catalyst 
-@parameters k₀₁ k₁₀ λ γ
 bursting_rn = @reaction_network begin 
     k₀₁, G0 --> G1 
     k₁₀, G1 --> G0 
     λ, G1 --> G1 + mRNA 
     γ, mRNA --> ∅
-end k₀₁ k₁₀ λ γ
+end 
 
 parameter_values = [k₀₁ => 0.05, k₁₀ => 0.1, λ => 5.0, γ => 0.5]
 model_from_catalyst = CmeModel(bursting_rn, parameter_values)
@@ -69,3 +68,4 @@ end
 
 
 
+
diff --git a/examples/toggleswitch_cme.jl b/examples/toggleswitch_fsp_variants.jl
similarity index 100%
rename from examples/toggleswitch_cme.jl
rename to examples/toggleswitch_fsp_variants.jl
diff --git a/src/cmemodel/catalyst_interface.jl b/src/cmemodel/catalyst_interface.jl
index e8ffcd1..98f5d16 100644
--- a/src/cmemodel/catalyst_interface.jl
+++ b/src/cmemodel/catalyst_interface.jl
@@ -5,9 +5,9 @@ function CmeModel(rn::Catalyst.ReactionSystem, parameter_values)
 
     species_ids = Catalyst.speciesmap(rn)
     parameter_ids = Catalyst.paramsmap(rn)
-    𝔖 = Catalyst.netstoichmat(rn)
+    stoich_matrix = Catalyst.netstoichmat(rn)
 
-    jump_rate_laws = [Catalyst.jumpratelaw(eq) for eq in rn.eqs]
+    jump_rate_laws = [Catalyst.jumpratelaw(eq) for eq in Catalyst.get_eqs(rn)]
 
     ## Convert Catalyst jump rate laws into 
     @variables t x[1:species_count] p[1:parameter_count]
@@ -37,5 +37,5 @@ function CmeModel(rn::Catalyst.ReactionSystem, parameter_values)
         end
     end
 
-    return CmeModel(𝔖, converted_propensities, parvec)
+    return CmeModel(stoich_matrix, converted_propensities, parvec)
 end
\ No newline at end of file
diff --git a/src/forwardsenscme/sparse/forwardsenscmesparse.jl b/src/forwardsenscme/sparse/forwardsenscmesparse.jl
index d242d97..b896d5a 100644
--- a/src/forwardsenscme/sparse/forwardsenscmesparse.jl
+++ b/src/forwardsenscme/sparse/forwardsenscmesparse.jl
@@ -145,7 +145,9 @@ function solve(model::CmeModelWithSensitivity,
         function affect!(integrator)
             DE.terminate!(integrator)
         end
-        # Error constraint for the intermediate solutions. The intermediate sinks must not grow beyond a linear function of time that reaches `fsptol` at the end time
+        # Error constraint for the intermediate solutions. 
+        # The intermediate sinks must not grow beyond a linear function 
+        # of time that reaches `fsptol` at the end time
         function fsp_error_constraint(u, t, integrator)
             sinks = u[n-sink_count+1:n]
             return sum(sinks) - fsptol * t / tend
@@ -159,8 +161,8 @@ function solve(model::CmeModelWithSensitivity,
             abstol = eps()
         )
 
-        fsensfspprob = DE.ODEProblem(sensfsprhs!, unow, (tnow, tend), p = get_parameters(model))
-        integrator = DE.init(fsensfspprob, sensfspalgorithm.ode_method, atol = odeatol, rtol = odertol, callback = fsp_cb, saveat = saveat, progress=true)
+        fsensfspprob = DE.ODEProblem(sensfsprhs!, unow, (tnow, tend), get_parameters(model))
+        integrator = DE.init(fsensfspprob, sensfspalgorithm.ode_method, abstol = odeatol, reltol = odertol, callback = fsp_cb, saveat = saveat, progress=true)
         DE.step!(integrator, tend - tnow, true)        
 
         for (t, u) in zip(integrator.sol.t, integrator.sol.u)
diff --git a/src/fspmatrix/fspmatrix.jl b/src/fspmatrix/fspmatrix.jl
index c7d72a1..9216add 100644
--- a/src/fspmatrix/fspmatrix.jl
+++ b/src/fspmatrix/fspmatrix.jl
@@ -1,6 +1,11 @@
-export FspMatrix, matvec, matvec!
+export AbstractFspMatrix, matvec, matvec!
+
+"""
+Abstract type for the FSP-truncated infinitesimal generator of the Chemical Master Equation.
+"""
+abstract type AbstractFspMatrix end
+
 
-abstract type FspMatrix end
 
 include("sparse/fspsparsematrix.jl")
 
diff --git a/src/fspmatrix/sparse/fspsparsematrix.jl b/src/fspmatrix/sparse/fspsparsematrix.jl
index 71b8b1c..4423522 100644
--- a/src/fspmatrix/sparse/fspsparsematrix.jl
+++ b/src/fspmatrix/sparse/fspsparsematrix.jl
@@ -6,7 +6,7 @@ export FspMatrixSparse, get_parameters, get_states, get_rowcount, get_colcount,
 """
 The transition rate matrix of the finite Markov chain approximation to the CME (including sink states). This type is based on the sparse representation of the truncated state space.
 """
-Base.@kwdef mutable struct FspMatrixSparse{NS,NR,IntT<:Integer,RealT<:AbstractFloat} <: FspMatrix
+Base.@kwdef mutable struct FspMatrixSparse{NS,NR,IntT<:Integer,RealT<:AbstractFloat} <: AbstractFspMatrix
     parameters::Vector{Any}
     states::Vector{MVector{NS,IntT}}
 
@@ -185,7 +185,9 @@ function size(A::FspMatrixSparse, dim::Integer)
     return (dim == 1) ? A.rowcount : A.colcount
 end
 
-# Matrix-vector multiplications
+#= 
+Matrix-vector multiplications
+=#
 """
 matvec!(out::Vector{RealT}, t::AbstractFloat, A::FspMatrixSparse, v::Vector{RealT})
 
@@ -195,12 +197,10 @@ function matvec!(out, t, A::FspMatrixSparse, v)
     if A.timeinvariant_matrix isa Nothing
         out .= 0.0
     else        
-        # out .= A.timeinvariant_matrix*v
         SArrays.mul!(out, A.timeinvariant_matrix, v)
     end
     θ = A.parameters
     for (i, r) in enumerate(A.separabletv_propensity_ids)
-        # out .+= A.propensities[r].tfactor(t, θ...)*A.separabletv_factormatrices[i]*v
         SArrays.mul!(out, A.separabletv_factormatrices[i], v, A.propensities[r].tfactor(t, θ), 1)        
     end    
     needupdate = false
@@ -211,7 +211,6 @@ function matvec!(out, t, A::FspMatrixSparse, v)
     for (i, r) in enumerate(A.jointtv_propensity_ids)
         (needupdate) && _update_sparsematrix!(A.jointtv_matrices[i], A.states, A.propensities[r],
             t, θ)
-        # out .+= A.jointtv_matrices[i]*v 
         SArrays.mul!(out, A.jointtv_matrices[i], v, 1.0, 1.0)                        
     end
     return nothing
@@ -258,18 +257,14 @@ function matvec(t, A::FspMatrixSparse, v)
     return w
 end
 
-function (A::FspMatrixSparse)(t::AbstractFloat)
-    return (t, A)
-end
-
 import Base: *
-function *(A_at_t::Tuple{AbstractFloat,FspMatrixSparse}, v::Vector{RealT}) where {RealT<:AbstractFloat}
-    t = A_at_t[1]
-    A = A_at_t[2]
-    return matvec(t, A, v)
-end
 
 function *(A::FspMatrixSparse, v::Vector{RealT}) where {RealT<:AbstractFloat}
     return matvec(0.0, A, v)
 end
 
+#=
+Interface to DifferentialEquations.jl's Linear operator 
+=#
+# mutable struct FspDiffEqOperator <: AbstractDiffEqOperator
+# end
diff --git a/src/fspvector/fspvector.jl b/src/fspvector/fspvector.jl
index 961892b..1198211 100644
--- a/src/fspvector/fspvector.jl
+++ b/src/fspvector/fspvector.jl
@@ -21,7 +21,7 @@ get_statedict(v::FspVectorSparse) = v.state2idx
 nnz(p::FspVectorSparse) = length(get_states(p))
 
 """
-    FspVectorSparse(states::Vector{MVector{NS, IntT}}, values::Vector{RealT})
+    FspVectorSparse(states::Vector{MVector{NS,IntT}}, values::Vector{RealT}; checksizes::Bool = true) where {NS,IntT<:Integer,RealT<:AbstractFloat}
 
 Construct a sparse representation of a FSP-truncated probability distribution with support `states` and proability values `values`.
 """
@@ -35,7 +35,7 @@ function FspVectorSparse(states::Vector{MVector{NS,IntT}}, values::Vector{RealT}
 end
 
 """
-    FspVectorSparse(statespace::AbstractStateSpaceSparse{NS, NR, IntT}, statevalpairs::Vector{Pair{VecT, RealT}})
+    FspVectorSparse(statespace::AbstractStateSpaceSparse{NS,NR,IntT}, statevalpairs::Vector{Pair{VecT,RealT}}) where {NS,NR,IntT<:Integer,VecT<:AbstractVector,RealT<:AbstractFloat}
 
 Construct a sparse representation of a FSP-truncated probability distribution on the state space `statespace` and proability values `values`.
 """
diff --git a/src/transientcme/sparse/rstepadapters.jl b/src/transientcme/sparse/rstepadapters.jl
index b269c60..c7c2561 100644
--- a/src/transientcme/sparse/rstepadapters.jl
+++ b/src/transientcme/sparse/rstepadapters.jl
@@ -1,3 +1,4 @@
+using SciMLBase: get_du!
 
 export RStepAdapter, SelectiveRStepAdapter
 
@@ -97,9 +98,9 @@ function adapt!(statespace::StateSpaceSparse, adapter::SelectiveRStepAdapter,
         deleteat!(p, dropids)
     end
     sink_count = get_sink_count(statespace)
-    du = similar(integrator.u)
-    f! = integrator.f
-    f!(du, integrator.u, [], t)
+    du = similar(integrator.u)    
+    get_du!(du, integrator)
+
     dsinks = du[end-sink_count+1:end]
     expandreactions = findall(dsinks .> 0)
     nold = get_state_count(statespace)
diff --git a/test/test_catalyst_interface.jl b/test/test_catalyst_interface.jl
index 188e68d..c119b32 100644
--- a/test/test_catalyst_interface.jl
+++ b/test/test_catalyst_interface.jl
@@ -27,13 +27,13 @@ end
 
 cmemodel1 = CmeModel(𝕊, [a1,a2,a3,a4], θ)
 
-@parameters k01 k10 α γ
+# @parameters k01 k10 α γ
 rn = @reaction_network begin
     k01, G0 --> G1
     k10, G1 --> G0
     α, G1 --> G1 + RNA
     γ*max(0.0, 1.0 - sin(π * t / L)), RNA --> ∅
-end k01 k10 α γ L
+end 
 
 cmemodel2 = CmeModel(rn, θ)
 @test get_species_count(cmemodel2) == 3
@@ -50,7 +50,7 @@ expand!(test_space, 100)
 pass = true 
 for t in test_times 
     for state in get_states(test_space)
-        pass &= cmemodel1.propensities[1](state,θ) ≈ cmemodel2.propensities[1](state,θ)
+        global pass &= cmemodel1.propensities[1](state,θ) ≈ cmemodel2.propensities[1](state,θ)
         pass &= cmemodel1.propensities[2](state,θ) ≈ cmemodel2.propensities[2](state,θ)
         pass &= cmemodel1.propensities[3](state,θ) ≈ cmemodel2.propensities[3](state,θ)
         pass &= cmemodel1.propensities[4](t,state,θ) ≈ cmemodel2.propensities[4](t,state,θ)
diff --git a/test/test_exponential.jl b/test/test_exponential.jl
new file mode 100644
index 0000000..6898e61
--- /dev/null
+++ b/test/test_exponential.jl
@@ -0,0 +1,26 @@
+# This test suite verifies that DifferentialEquations.jl Exponential integrators can be used with the FSP 
+using NumCME
+using DifferentialEquations
+using StaticArrays 
+import Catalyst 
+
+## Test Fixture: telegraph gene expression
+Catalyst.@parameters k₀₁ k₁₀ λ γ
+bursting_rn = Catalyst.@reaction_network begin 
+    k₀₁, G0 --> G1 
+    k₁₀, G1 --> G0 
+    λ, G1 --> G1 + mRNA 
+    γ, mRNA --> ∅
+end k₀₁ k₁₀ λ γ
+
+parameter_values = [k₀₁ => 0.05, k₁₀ => 0.1, λ => 5.0, γ => 0.5]
+model_from_catalyst = CmeModel(bursting_rn, parameter_values)
+
+fspalgorithm = AdaptiveFspSparse(
+    ode_method=LinearExponential(),
+    space_adapter=RStepAdapter(5, 10, true)
+)
+p0 = FspVectorSparse([@MVector [1, 0, 0]], [1.0])
+fspsol = solve(model_from_catalyst, p0, (0.0, 100.0), fspalgorithm; saveat=[0.0, 50.0, 100.0])
+
+
diff --git a/test/test_fspmat.jl b/test/test_fspmat.jl
index d5301b3..71feb07 100644
--- a/test/test_fspmat.jl
+++ b/test/test_fspmat.jl
@@ -42,9 +42,9 @@ expand!(𝔛, 2)
 @test size(𝐀, 1) == get_state_count(𝔛) + get_sink_count(𝔛)
 @test size(𝐀, 2) == get_state_count(𝔛) + get_sink_count(𝔛)
 𝐯 = ones(Float64, size(𝐀, 1))
-𝐰 = 𝐀(1.0) * 𝐯
+𝐰 = matvec(1.0, 𝐀, 𝐯)
 @test sum(𝐰) ≈ 0.0 atol = 1.0e-14
-𝐰 = 𝐀 * 𝐯
+𝐰 = matvec(0.0, 𝐀, 𝐯)
 @test sum(𝐰) ≈ 0.0 atol = 1.0e-14
 
 # Test mat-vec for time-varying matrix
@@ -54,15 +54,15 @@ A1 = FspMatrixSparse(𝔛, propensities_tv, parameters=θ)
 @test size(A1, 1) == get_state_count(𝔛) + get_sink_count(𝔛)
 @test size(A1, 2) == get_state_count(𝔛) + get_sink_count(𝔛)
 𝐯 = ones(Float64, size(A1, 1))
-w1 = A1(1.0) * 𝐯
+w1 = matvec(1.0, A1, 𝐯)
 @test sum(w1) ≈ 0.0 atol = 1.0e-14
-w1 = A1 * 𝐯
+w1 = matvec(0.0, A1, 𝐯)
 @test sum(w1) ≈ 0.0 atol = 1.0e-14
 
 A2 = FspMatrixSparse(𝔛, propensities_tvj, parameters=θ)
-w2 = A2(1.0) * 𝐯
+w2 = matvec(1.0, A2, 𝐯)
 @test sum(w2) ≈ 0.0 atol = 1.0e-14
-w2 = A2 * 𝐯
+w2 = matvec(0.0, A2, 𝐯)
 
 @test sum(w2) ≈ 0.0 atol = 1.0e-14
 @test norm(w1 -w2) ≈ 0