diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 6b48ad2..0cddec6 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -1,41 +1,40 @@
-name: CI
-on:
-  - push
-  - pull_request
-jobs:
-  test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
-    runs-on: ${{ matrix.os }}
-    strategy:
-      fail-fast: false
-      matrix:
-        version:
-          - '1.5'
-          - '1.6'
-          - '1.7'
-        os:
-          - ubuntu-latest
-        arch:
-          - x64
-    steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
-        with:
-          version: ${{ matrix.version }}
-          arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
-      - uses: julia-actions/julia-buildpkg@v1
-      - uses: julia-actions/julia-runtest@v1
-      - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
-        with:
-          file: lcov.info
+name: CI
+on:
+  - push
+  - pull_request
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.6'
+          - '1.9'
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v1
+        with:
+          file: lcov.info
diff --git a/Project.toml b/Project.toml
index 720c8fa..bb42e86 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "FiniteStateProjection"
 uuid = "069e79ea-d681-44e8-b935-95bdaf9e8f28"
 authors = ["kaandocal"]
-version = "0.2.1"
+version = "0.3.0"
 
 [deps]
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
@@ -12,12 +12,12 @@ RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-Catalyst = "9, 10, 11, 12"
+Catalyst = "13"
 DiffEqBase = "6"
 MacroTools = "^0.5.5"
 Reexport = "1"
 RuntimeGeneratedFunctions = "0.5"
-julia = "1.5"
+julia = "1.6"
 
 [extras]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
diff --git a/README.md b/README.md
index adc9845..5e797ba 100644
--- a/README.md
+++ b/README.md
@@ -9,7 +9,7 @@ Finite State Projection [[1]](#1)  algorithms for chemical reaction networks bas
 - FSP equations are generated as `ODEFunction`/`ODEProblem`s and can be solved with [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl), with on-the-fly generation of targeted functions for improved performance
 - The Chemical Master Equation can be represented as a `SparseMatrixCSC`
 
-More information is available in the [documentation](https://kaandocal.github.io/FiniteStateProjection.jl/dev/). Please feel free to open issues and submit pull requests! 
+More information is available in the [documentation](https://kaandocal.github.io/FiniteStateProjection.jl/dev/). Please feel free to open issues and submit pull requests!
 
 ## Examples
 ### Birth-Death System
@@ -20,7 +20,7 @@ using OrdinaryDiffEq
 rn = @reaction_network begin
     σ, 0 --> A
     d, A --> 0
-end σ d
+end
 
 sys = FSPSystem(rn)
 
@@ -30,7 +30,7 @@ ps = [ 10.0, 1.0 ]
 # Initial distribution (over 1 species)
 # Here we start with 0 copies of A
 u0 = zeros(50)
-u0[1] = 1.0 
+u0[1] = 1.0
 
 prob = convert(ODEProblem, sys, u0, (0, 10.0), ps)
 sol = solve(prob, Vern7())
@@ -47,7 +47,7 @@ rn = @reaction_network begin
     σ_off, G_on --> 0
     ρ, G_on --> G_on + M
     d, M --> 0
-end σ_on σ_off ρ d
+end
 
 sys = FSPSystem(rn)
 
diff --git a/demo/birth_death.jl b/demo/birth_death.jl
index e5c5c7c..9f3247d 100644
--- a/demo/birth_death.jl
+++ b/demo/birth_death.jl
@@ -7,7 +7,7 @@ using PyPlot
 rs = @reaction_network begin
     r1, 0 --> A
     r2, A --> 0
-end r1 r2
+end
 
 ##
 
diff --git a/demo/telegraph.jl b/demo/telegraph.jl
index a42c3d2..1bcc820 100644
--- a/demo/telegraph.jl
+++ b/demo/telegraph.jl
@@ -9,7 +9,7 @@ rs = @reaction_network begin
     r2, G_on --> 0
     r3, G_on --> G_on + M
     r4, M --> 0
-end r1 r2 r3 r4
+end
 
 ##
 
@@ -20,9 +20,9 @@ ps = [ 0.25, 0.15, 15.0, 1.0 ]
 
 # Initial values
 # Since G_on + G_off = const. we do not have to model the two
-# separately. Use reduced_species(sys) to get the list of 
+# separately. Use reduced_species(sys) to get the list of
 # species we actually have to model:
-# 
+#
 # julia> reduced_species(sys)
 # 2-element Vector{Int64}:
 #  1
diff --git a/docs/Project.toml b/docs/Project.toml
index dfa65cd..3a52a5d 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,2 +1,5 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+
+[compat]
+Documenter = "0.27"
diff --git a/docs/src/examples.md b/docs/src/examples.md
index 09dbee4..45f7655 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -11,7 +11,7 @@ using OrdinaryDiffEq
 rn = @reaction_network begin
     σ, 0 --> A
     d, A --> 0
-end σ d
+end
 
 sys = FSPSystem(rn)
 
@@ -21,7 +21,7 @@ ps = [ 10.0, 1.0 ]
 # Initial distribution (over 1 species)
 # Here we start with 0 copies of A
 u0 = zeros(50)
-u0[1] = 1.0 
+u0[1] = 1.0
 
 prob = ODEProblem(sys, u0, (0, 10.0), ps)
 sol = solve(prob, Vern7())
@@ -47,7 +47,7 @@ rn = @reaction_network begin
     σ_off, G_on --> 0
     ρ, G_on --> G_on + M
     d, M --> 0
-end σ_on σ_off ρ d
+end
 
 sys = FSPSystem(rn)
 
diff --git a/docs/src/troubleshoot.md b/docs/src/troubleshoot.md
index 3f92b36..6e69589 100644
--- a/docs/src/troubleshoot.md
+++ b/docs/src/troubleshoot.md
@@ -42,7 +42,7 @@ This point might seem obvious, but errors in the rate functions, or an incorrect
 rn = @reaction_network begin
    σ * (N - I), I --> 2I
    ρ, I --> 0
-end σ ρ N
+end
 
 sys_fsp = FSPSystem(rn)
 ```
diff --git a/src/fspsystem.jl b/src/fspsystem.jl
index bc45153..a2081e7 100644
--- a/src/fspsystem.jl
+++ b/src/fspsystem.jl
@@ -1,5 +1,5 @@
-""" 
-Thin wrapper around `Catalyst.ReactionSystem` for use with this package. 
+"""
+Thin wrapper around `Catalyst.ReactionSystem` for use with this package.
 
 Constructor: `FSPSystem(rs::ReactionSystem[, ih=DefaultIndexHandler(); combinatoric_ratelaw::Bool=true])`
 """
@@ -10,22 +10,27 @@ struct FSPSystem{IHT <: AbstractIndexHandler, RT}
 end
 
 function FSPSystem(rs::ReactionSystem, ih=DefaultIndexHandler(); combinatoric_ratelaw::Bool=true)
+    isempty(Catalyst.get_systems(rs)) ||
+        error("Supported Catalyst models can not contain subsystems. Please use `rs = Catalyst.flatten(rs::ReactionSystem)` to generate a single system with no subsystems from your Catalyst model.")
+    any(eq -> !(eq isa Reaction), equations(rs)) &&
+        error("Catalyst models that include constraint ODEs or algebraic equations are not supported.")
+
     rfs = create_ratefuncs(rs, ih; combinatoric_ratelaw=combinatoric_ratelaw)
     FSPSystem(rs, ih, rfs)
 end
 
-""" 
+"""
     build_ratefuncs(rs, ih; state_sym::Symbol, combinatoric_ratelaw::Bool)::Vector
 
-Return the rate functions converted to Julia expressions in the state variable 
+Return the rate functions converted to Julia expressions in the state variable
 `state_sym`. Abundances of the species are computed using `getsubstitutions`.
 
 See also: [`getsubstitutions`](@ref), [`build_rhs`](@ref)
 """
 function build_ratefuncs(rs::ReactionSystem, ih::AbstractIndexHandler; state_sym::Symbol, combinatoric_ratelaw::Bool=true)
     substitutions = getsubstitutions(ih, rs, state_sym=state_sym)
-    
-    return map(Catalyst.get_eqs(rs)) do reac
+
+    return map(Catalyst.reactions(rs)) do reac
         jrl = jumpratelaw(reac; combinatoric_ratelaw)
         jrl_s = substitute(jrl, substitutions)
         toexpr(jrl_s)
@@ -34,12 +39,12 @@ end
 
 function create_ratefuncs(rs::ReactionSystem, ih::AbstractIndexHandler; combinatoric_ratelaw::Bool=true)
     paramsyms = Symbol.(Catalyst.parameters(rs))
-    
-    return tuple(map(ex -> compile_ratefunc(ex, paramsyms), 
+
+    return tuple(map(ex -> compile_ratefunc(ex, paramsyms),
                      build_ratefuncs(rs, ih; state_sym=:idx_in, combinatoric_ratelaw))...)
-end 
+end
 
-function compile_ratefunc(ex_rf, params) 
+function compile_ratefunc(ex_rf, params)
     # Make this nicer in the future
     ex = :((idx_in, t, $(params...)) -> $(ex_rf)) |> MacroTools.flatten
     @RuntimeGeneratedFunction(ex)
diff --git a/src/indexhandlers.jl b/src/indexhandlers.jl
index b25ce7d..e310ae6 100644
--- a/src/indexhandlers.jl
+++ b/src/indexhandlers.jl
@@ -112,5 +112,7 @@ end
 Defines the abundance of species ``S_i`` to be `state_sym[i] - offset`.
 """
 function getsubstitutions(ih::DefaultIndexHandler, rs::ReactionSystem; state_sym::Symbol)
-    Dict(symbol => Term{Number}(Base.getindex, (state_sym, i)) - ih.offset for (i, symbol) in enumerate(species(rs)))
+    nspecs = numspecies(rs)
+    state_sym_vec = ModelingToolkit.value.(ModelingToolkit.scalarize((@variables ($state_sym)[1:nspecs])[1]))
+    Dict(symbol => state_sym_vec[i] - ih.offset for (i, symbol) in enumerate(species(rs)))
 end
diff --git a/src/matrix.jl b/src/matrix.jl
index 588ba02..ddd19c0 100644
--- a/src/matrix.jl
+++ b/src/matrix.jl
@@ -6,7 +6,7 @@ function create_sparsematrix(sys::FSPSystem, dims::NTuple, ps, t)
     J = Int[]
     V = Float64[]
 
-    predsize = Ntot * (length(Catalyst.get_eqs(sys.rs)) + 1)
+    predsize = Ntot * (length(Catalyst.equations(sys.rs)) + 1)
 
     sizehint!(I, predsize)
     sizehint!(J, predsize)
@@ -33,7 +33,7 @@ function create_sparsematrix(sys::FSPSystem, dims::NTuple, ps, t)
             push!(I, lind[idx_cout])
             push!(J, lind[idx_cin])
 
-            rate = rf(idx_cin, t, ps...) 
+            rate = rf(idx_cin, t, ps...)
             push!(V, rate)
         end
     end
@@ -49,7 +49,7 @@ function create_sparsematrix_ss(sys::FSPSystem, dims::NTuple, ps)
     J = Int[]
     V = Float64[]
 
-    predsize = 2 * Ntot * length(Catalyst.get_eqs(sys.rs))
+    predsize = 2 * Ntot * length(Catalyst.equations(sys.rs))
 
     sizehint!(I, predsize)
     sizehint!(J, predsize)
@@ -79,8 +79,8 @@ end
 """
     Base.convert(::Type{SparseMatrixCSC}, sys::FSPSystem, dims::NTuple, ps, t::Real)
 
-Convert the reaction system into a sparse matrix defining the right-hand side of the 
-Chemical Master Equation. `dims` is a tuple denoting the dimensions of the FSP and 
+Convert the reaction system into a sparse matrix defining the right-hand side of the
+Chemical Master Equation. `dims` is a tuple denoting the dimensions of the FSP and
 `ps` is the tuple of parameters. The sparse matrix works on the flattened version
 of the state obtained using `vec`.
 """
@@ -91,7 +91,7 @@ end
 """
     Base.convert(::Type{SparseMatrixCSC}, sys::FSPSystem, dims::NTuple, ps, ::SteadyState)
 
-Convert the reaction system into a sparse matrix defining the right-hand side of the 
+Convert the reaction system into a sparse matrix defining the right-hand side of the
 Chemical Master Equation, steady-state version.
 """
 function Base.convert(::Type{SparseMatrixCSC}, sys::FSPSystem, dims::NTuple, ps, ::SteadyState)
diff --git a/test/birthdeath2D.jl b/test/birthdeath2D.jl
index 03c25ac..ce92d4c 100644
--- a/test/birthdeath2D.jl
+++ b/test/birthdeath2D.jl
@@ -9,28 +9,27 @@ using Sundials
 
 marg(vec; dims) = dropdims(sum(vec; dims); dims)
 
-@parameters r1, r2, s1, s2
 rs = @reaction_network begin
     r1, 0 --> A
     r2, A --> 0
     s1, 0 --> B
     s2, B --> 0
-end r1 r2 s1 s2
+end
 
 sys = FSPSystem(rs)
 
 prs = exp.(2 .* rand(2))
-ps = [ prs[1], prs[1] / exp(4 * rand()), 
+ps = [ prs[1], prs[1] / exp(4 * rand()),
        prs[2], prs[2] / exp(4 * rand()) ]
 
-Nmax = 130
+Nmax = 100
 u0 = zeros(Nmax+1, Nmax+1)
 u0[1] = 1.0
 
 tt = [ 0.25, 1.0, 10.0 ]
 
 prob = convert(ODEProblem, sys, u0, 10.0, ps)
-sol = solve(prob, Vern7(), abstol=1e-9, reltol=1e-6, saveat=tt)
+sol = solve(prob, Vern7(), abstol=1e-6, saveat=tt)
 
 @test marg(sol.u[1], dims=2) ≈ pdf.(Poisson(ps[1] / ps[2] * (1 - exp(-ps[2] * tt[1]))), 0:Nmax) atol=1e-4
 @test marg(sol.u[1], dims=1) ≈ pdf.(Poisson(ps[3] / ps[4] * (1 - exp(-ps[4] * tt[1]))), 0:Nmax) atol=1e-4
@@ -45,7 +44,7 @@ A = convert(SparseMatrixCSC, sys, (Nmax+1, Nmax+1), ps, 0)
 f = (du,u,t) -> mul!(du, A, u)
 
 probA = ODEProblem(f, u0, 10.0)
-solA = solve(prob, Vern7(), abstol=1e-9, reltol=1e-6, saveat=tt)
+solA = solve(prob, Vern7(), abstol=1e-6, saveat=tt)
 
 @test sol.u[1] ≈ solA.u[1] atol=1e-4
 @test sol.u[2] ≈ solA.u[2] atol=1e-4
diff --git a/test/feedbackloop.jl b/test/feedbackloop.jl
index b9d3976..973e7ac 100644
--- a/test/feedbackloop.jl
+++ b/test/feedbackloop.jl
@@ -6,15 +6,13 @@ using SteadyStateDiffEq
 using FiniteStateProjection
 using SparseArrays
 
-@parameters σ_off σ_on ρ_on ρ_off d
-
 rs = @reaction_network begin
 	σ_off, G + P → 0
 	σ_on * (1 - G), 0 ⇒ G + P
 	ρ_on, G → G + P
 	ρ_off * (1-G), 0 ⇒ P
 	d, P → 0
-end σ_off σ_on ρ_on ρ_off d
+end
 
 ps = [ 1, 0.1, 20, 1, 1]
 Nmax = 200
@@ -49,12 +47,12 @@ u0 = zeros(2, Nmax)
 u0[1] = 1.0
 
 prob = convert(ODEProblem, sys, u0, maximum(tt), ps)
-sol = solve(prob, Vern7(), abstol=1e-9, reltol=1e-6, saveat=tt)
+sol = solve(prob, Vern7(), abstol=1e-6, saveat=tt)
 
 f = (du,u,t) -> mul!(du, A, u)
 
 probA = ODEProblem(f, u0, 10.0)
-solA = solve(prob, Vern7(), abstol=1e-9, reltol=1e-6, saveat=tt)
+solA = solve(prob, Vern7(), abstol=1e-6, saveat=tt)
 
 @test sol.u[1] ≈ solA.u[1] atol=1e-4
 @test sol.u[2] ≈ solA.u[2] atol=1e-4