diff --git a/Project.toml b/Project.toml
index b2296d4f..998f2e62 100644
--- a/Project.toml
+++ b/Project.toml
@@ -39,7 +39,7 @@ RecipesBase = "0.7, 0.8, 1.0"
 Requires = "1.0"
 StaticArraysCore = "1.1"
 Statistics = "1"
-SymbolicIndexingInterface = "0.1, 0.2"
+SymbolicIndexingInterface = "0.3"
 Tables = "1"
 Zygote = "0.6.56"
 julia = "1.6"
diff --git a/docs/make.jl b/docs/make.jl
index eb6f13df..221c6057 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -6,13 +6,13 @@ cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true)
 include("pages.jl")
 
 makedocs(sitename = "RecursiveArrayTools.jl",
-         authors = "Chris Rackauckas",
-         modules = [RecursiveArrayTools],
-         clean = true, doctest = false, linkcheck = true,
-         warnonly = [:missing_docs],
-         format = Documenter.HTML(assets = ["assets/favicon.ico"],
-                                  canonical = "https://docs.sciml.ai/RecursiveArrayTools/stable/"),
-         pages = pages)
+    authors = "Chris Rackauckas",
+    modules = [RecursiveArrayTools],
+    clean = true, doctest = false, linkcheck = true,
+    warnonly = [:missing_docs],
+    format = Documenter.HTML(assets = ["assets/favicon.ico"],
+        canonical = "https://docs.sciml.ai/RecursiveArrayTools/stable/"),
+    pages = pages)
 
 deploydocs(repo = "github.com/SciML/RecursiveArrayTools.jl.git";
-           push_preview = true)
+    push_preview = true)
diff --git a/ext/RecursiveArrayToolsMeasurementsExt.jl b/ext/RecursiveArrayToolsMeasurementsExt.jl
index 50eaaa24..d631655e 100644
--- a/ext/RecursiveArrayToolsMeasurementsExt.jl
+++ b/ext/RecursiveArrayToolsMeasurementsExt.jl
@@ -4,8 +4,8 @@ import RecursiveArrayTools
 isdefined(Base, :get_extension) ? (import Measurements) : (import ..Measurements)
 
 function RecursiveArrayTools.recursive_unitless_bottom_eltype(a::Type{
-                                                                      <:Measurements.Measurement
-                                                                      })
+    <:Measurements.Measurement,
+})
     typeof(oneunit(a))
 end
 
diff --git a/ext/RecursiveArrayToolsMonteCarloMeasurementsExt.jl b/ext/RecursiveArrayToolsMonteCarloMeasurementsExt.jl
index 76825d6d..35250da6 100644
--- a/ext/RecursiveArrayToolsMonteCarloMeasurementsExt.jl
+++ b/ext/RecursiveArrayToolsMonteCarloMeasurementsExt.jl
@@ -1,15 +1,18 @@
 module RecursiveArrayToolsMonteCarloMeasurementsExt
 
 import RecursiveArrayTools
-isdefined(Base, :get_extension) ? (import MonteCarloMeasurements) : (import ..MonteCarloMeasurements)
+isdefined(Base, :get_extension) ? (import MonteCarloMeasurements) :
+(import ..MonteCarloMeasurements)
 
 function RecursiveArrayTools.recursive_unitless_bottom_eltype(a::Type{
-                                                                      <:MonteCarloMeasurements.Particles
-                                                                      })
+    <:MonteCarloMeasurements.Particles,
+})
     typeof(one(a))
 end
 
-function RecursiveArrayTools.recursive_unitless_eltype(a::Type{<:MonteCarloMeasurements.Particles})
+function RecursiveArrayTools.recursive_unitless_eltype(a::Type{
+    <:MonteCarloMeasurements.Particles,
+})
     typeof(one(a))
 end
 
diff --git a/ext/RecursiveArrayToolsTrackerExt.jl b/ext/RecursiveArrayToolsTrackerExt.jl
index a06163c2..845aff72 100644
--- a/ext/RecursiveArrayToolsTrackerExt.jl
+++ b/ext/RecursiveArrayToolsTrackerExt.jl
@@ -4,12 +4,12 @@ import RecursiveArrayTools
 isdefined(Base, :get_extension) ? (import Tracker) : (import ..Tracker)
 
 function RecursiveArrayTools.recursivecopy!(b::AbstractArray{T, N},
-                                            a::AbstractArray{T2, N}) where {
-                                                                            T <:
-                                                                            Tracker.TrackedArray,
-                                                                            T2 <:
-                                                                            Tracker.TrackedArray,
-                                                                            N}
+    a::AbstractArray{T2, N}) where {
+    T <:
+    Tracker.TrackedArray,
+    T2 <:
+    Tracker.TrackedArray,
+    N}
     @inbounds for i in eachindex(a)
         b[i] = copy(a[i])
     end
diff --git a/ext/RecursiveArrayToolsZygoteExt.jl b/ext/RecursiveArrayToolsZygoteExt.jl
index 644477fa..5f93bf96 100644
--- a/ext/RecursiveArrayToolsZygoteExt.jl
+++ b/ext/RecursiveArrayToolsZygoteExt.jl
@@ -14,7 +14,7 @@ end
 ChainRulesCore.ProjectTo(x::VectorOfArray) = ChainRulesCore.ProjectTo{VectorOfArray}()
 
 function ChainRulesCore.rrule(T::Type{<:RecursiveArrayTools.GPUArraysCore.AbstractGPUArray},
-                              xs::AbstractVectorOfArray)
+    xs::AbstractVectorOfArray)
     T(xs), ȳ -> (ChainRulesCore.NoTangent(), ȳ)
 end
 
@@ -28,7 +28,7 @@ end
 end
 
 @adjoint function getindex(VA::AbstractVectorOfArray,
-                           i::Union{BitArray, AbstractArray{Bool}})
+    i::Union{BitArray, AbstractArray{Bool}})
     function AbstractVectorOfArray_getindex_adjoint(Δ)
         Δ′ = [(i[j] ? Δ[j] : FillArrays.Fill(zero(eltype(x)), size(x)))
               for (x, j) in zip(VA.u, 1:length(VA))]
@@ -48,7 +48,7 @@ end
 end
 
 @adjoint function getindex(VA::AbstractVectorOfArray,
-                           i::Union{Int, AbstractArray{Int}})
+    i::Union{Int, AbstractArray{Int}})
     function AbstractVectorOfArray_getindex_adjoint(Δ)
         Δ′ = [(i[j] ? Δ[j] : FillArrays.Fill(zero(eltype(x)), size(x)))
               for (x, j) in zip(VA.u, 1:length(VA))]
@@ -65,8 +65,8 @@ end
 end
 
 @adjoint function getindex(VA::AbstractVectorOfArray, i::Int,
-                           j::Union{Int, AbstractArray{Int}, CartesianIndex,
-                                    Colon, BitArray, AbstractArray{Bool}}...)
+    j::Union{Int, AbstractArray{Int}, CartesianIndex,
+        Colon, BitArray, AbstractArray{Bool}}...)
     function AbstractVectorOfArray_getindex_adjoint(Δ)
         Δ′ = VectorOfArray([zero(x) for (x, j) in zip(VA.u, 1:length(VA))])
         Δ′[i, j...] = Δ
@@ -76,11 +76,11 @@ end
 end
 
 @adjoint function ArrayPartition(x::S,
-                                 ::Type{Val{copy_x}} = Val{false}) where {
-                                                                          S <:
-                                                                          Tuple,
-                                                                          copy_x
-                                                                          }
+    ::Type{Val{copy_x}} = Val{false}) where {
+    S <:
+    Tuple,
+    copy_x,
+}
     function ArrayPartition_adjoint(_y)
         y = Array(_y)
         starts = vcat(0, cumsum(reduce(vcat, length.(x))))
@@ -93,14 +93,21 @@ end
 
 @adjoint function VectorOfArray(u)
     VectorOfArray(u),
-    y -> (VectorOfArray([y[ntuple(x -> Colon(), ndims(y) - 1)..., i]
-                         for i in 1:size(y)[end]]),)
+    y -> begin
+        y isa Ref && (y = VectorOfArray(y[].u))
+        (VectorOfArray([y[ntuple(x -> Colon(), ndims(y.u) - 1)..., i]
+                        for i in 1:size(y.u)[end]]),)
+    end
 end
 
 @adjoint function DiffEqArray(u, t)
     DiffEqArray(u, t),
-    y -> (DiffEqArray([y[ntuple(x -> Colon(), ndims(y) - 1)..., i] for i in 1:size(y)[end]],
-                      t), nothing)
+    y -> begin
+        y isa Ref && (y = VectorOfArray(y[].u))
+        (DiffEqArray([y[ntuple(x -> Colon(), ndims(y.u) - 1)..., i]
+                       for i in 1:size(y.u)[end]],
+            t), nothing)
+    end
 end
 
 @adjoint function literal_getproperty(A::ArrayPartition, ::Val{:x})
diff --git a/src/RecursiveArrayTools.jl b/src/RecursiveArrayTools.jl
index 51fa3b0a..46e47fb0 100644
--- a/src/RecursiveArrayTools.jl
+++ b/src/RecursiveArrayTools.jl
@@ -6,14 +6,14 @@ module RecursiveArrayTools
 
 using DocStringExtensions
 using RecipesBase, StaticArraysCore, Statistics,
-      ArrayInterface, LinearAlgebra
+    ArrayInterface, LinearAlgebra
 using SymbolicIndexingInterface
 
 import Adapt
 
 import Tables, IteratorInterfaceExtensions
 
-abstract type AbstractVectorOfArray{T, N, A} <: AbstractArray{T, N} end
+abstract type AbstractVectorOfArray{T, N, A} end
 abstract type AbstractDiffEqArray{T, N, A} <: AbstractVectorOfArray{T, N, A} end
 
 include("utils.jl")
@@ -31,18 +31,24 @@ Base.convert(T::Type{<:GPUArraysCore.AbstractGPUArray}, VA::AbstractVectorOfArra
 import Requires
 @static if !isdefined(Base, :get_extension)
     function __init__()
-        Requires.@require Measurements="eff96d63-e80a-5855-80a2-b1b0885c5ab7" begin include("../ext/RecursiveArrayToolsMeasurementsExt.jl") end
-        Requires.@require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin include("../ext/RecursiveArrayToolsTrackerExt.jl") end
-        Requires.@require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" begin include("../ext/RecursiveArrayToolsZygoteExt.jl") end
+        Requires.@require Measurements="eff96d63-e80a-5855-80a2-b1b0885c5ab7" begin
+            include("../ext/RecursiveArrayToolsMeasurementsExt.jl")
+        end
+        Requires.@require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
+            include("../ext/RecursiveArrayToolsTrackerExt.jl")
+        end
+        Requires.@require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" begin
+            include("../ext/RecursiveArrayToolsZygoteExt.jl")
+        end
     end
 end
 
 export VectorOfArray, DiffEqArray, AbstractVectorOfArray, AbstractDiffEqArray,
-       AllObserved, vecarr_to_vectors, tuples
+    AllObserved, vecarr_to_vectors, tuples
 
 export recursivecopy, recursivecopy!, recursivefill!, vecvecapply, copyat_or_push!,
-       vecvec_to_mat, recursive_one, recursive_mean, recursive_bottom_eltype,
-       recursive_unitless_bottom_eltype, recursive_unitless_eltype
+    vecvec_to_mat, recursive_one, recursive_mean, recursive_bottom_eltype,
+    recursive_unitless_bottom_eltype, recursive_unitless_eltype
 
 export ArrayPartition
 
diff --git a/src/array_partition.jl b/src/array_partition.jl
index 62cd9948..2bfc777e 100644
--- a/src/array_partition.jl
+++ b/src/array_partition.jl
@@ -91,10 +91,14 @@ Base.zero(A::ArrayPartition, dims::NTuple{N, Int}) where {N} = zero(A)
 ## Array
 
 Base.Array(A::ArrayPartition) = reduce(vcat, Array.(A.x))
-function Base.Array(VA::AbstractVectorOfArray{T, N, A}) where {T, N,
-                                                               A <: AbstractVector{
-                                                                              <:ArrayPartition
-                                                                              }}
+function Base.Array(VA::AbstractVectorOfArray{
+    T,
+    N,
+    A,
+}) where {T, N,
+    A <: AbstractVector{
+        <:ArrayPartition,
+    }}
     reduce(hcat, Array.(VA.u))
 end
 
@@ -115,7 +119,7 @@ Base.ones(A::ArrayPartition, dims::NTuple{N, Int}) where {N} = ones(A)
 
 # mutable iff all components of ArrayPartition are mutable
 @generated function ArrayInterface.ismutable(::Type{<:ArrayPartition{T, S}}) where {T, S
-                                                                                    }
+}
     res = all(ArrayInterface.ismutable, S.parameters)
     return :($res)
 end
@@ -215,18 +219,6 @@ Base.@propagate_inbounds function Base.getindex(A::ArrayPartition, i::Int)
     end
 end
 
-"""
-    getindex(A::ArrayPartition, i::Int, j...)
-
-Returns the entry at index `j...` of the `i`th partition of `A`.
-"""
-Base.@propagate_inbounds function Base.getindex(A::ArrayPartition, i::Int, j...)
-    @boundscheck 0 < i <= length(A.x) || throw(BoundsError(A.x, i))
-    @inbounds b = A.x[i]
-    @boundscheck checkbounds(b, j...)
-    @inbounds return b[j...]
-end
-
 """
     getindex(A::ArrayPartition, i::Colon, j...)
 
@@ -256,7 +248,7 @@ end
 
 # workaround for https://github.com/SciML/RecursiveArrayTools.jl/issues/49
 function Base._unsafe_getindex(::IndexStyle, A::ArrayPartition,
-                               I::Vararg{Union{Real, AbstractArray}, N}) where {N}
+    I::Vararg{Union{Real, AbstractArray}, N}) where {N}
     # This is specifically not inlined to prevent excessive allocations in type unstable code
     shape = Base.index_shape(I...)
     dest = similar(A.x[1], shape)
@@ -264,18 +256,10 @@ function Base._unsafe_getindex(::IndexStyle, A::ArrayPartition,
     return dest
 end
 
-Base._maybe_reshape(::IndexCartesian, A::ArrayPartition, I::Vararg{Union{Real, AbstractArray}, N}) where {N} = Vector(A)
-
-"""
-    setindex!(A::ArrayPartition, v, i::Int, j...)
-
-Set the entry at index `j...` of the `i`th partition of `A` to `v`.
-"""
-Base.@propagate_inbounds function Base.setindex!(A::ArrayPartition, v, i::Int, j...)
-    @boundscheck 0 < i <= length(A.x) || throw(BoundsError(A.x, i))
-    @inbounds b = A.x[i]
-    @boundscheck checkbounds(b, j...)
-    @inbounds b[j...] = v
+function Base._maybe_reshape(::IndexCartesian,
+    A::ArrayPartition,
+    I::Vararg{Union{Real, AbstractArray}, N}) where {N}
+    Vector(A)
 end
 
 ## recursive methods
@@ -321,19 +305,19 @@ end
 
 # promotion rules
 @inline function Broadcast.BroadcastStyle(::ArrayPartitionStyle{AStyle},
-                                          ::ArrayPartitionStyle{BStyle}) where {AStyle,
-                                                                                BStyle}
+    ::ArrayPartitionStyle{BStyle}) where {AStyle,
+    BStyle}
     ArrayPartitionStyle(Broadcast.BroadcastStyle(AStyle(), BStyle()))
 end
 function Broadcast.BroadcastStyle(::ArrayPartitionStyle{Style},
-                                  ::Broadcast.DefaultArrayStyle{0}) where {
-                                                                           Style <:
-                                                                           Broadcast.BroadcastStyle
-                                                                           }
+    ::Broadcast.DefaultArrayStyle{0}) where {
+    Style <:
+    Broadcast.BroadcastStyle,
+}
     ArrayPartitionStyle{Style}()
 end
 function Broadcast.BroadcastStyle(::ArrayPartitionStyle,
-                                  ::Broadcast.DefaultArrayStyle{N}) where {N}
+    ::Broadcast.DefaultArrayStyle{N}) where {N}
     Broadcast.DefaultArrayStyle{N}()
 end
 
@@ -343,11 +327,11 @@ combine_styles(args::Tuple{}) = Broadcast.DefaultArrayStyle{0}()
 end
 @inline function combine_styles(args::Tuple{Any, Any})
     Broadcast.result_style(Broadcast.BroadcastStyle(args[1]),
-                           Broadcast.BroadcastStyle(args[2]))
+        Broadcast.BroadcastStyle(args[2]))
 end
 @inline function combine_styles(args::Tuple)
     Broadcast.result_style(Broadcast.BroadcastStyle(args[1]),
-                           combine_styles(Base.tail(args)))
+        combine_styles(Base.tail(args)))
 end
 
 function Broadcast.BroadcastStyle(::Type{ArrayPartition{T, S}}) where {T, S}
@@ -355,9 +339,11 @@ function Broadcast.BroadcastStyle(::Type{ArrayPartition{T, S}}) where {T, S}
     ArrayPartitionStyle(Style)
 end
 
-@inline function Base.copy(bc::Broadcast.Broadcasted{ArrayPartitionStyle{Style}}) where {
-                                                                                         Style
-                                                                                         }
+@inline function Base.copy(bc::Broadcast.Broadcasted{
+    ArrayPartitionStyle{Style},
+}) where {
+    Style,
+}
     N = npartitions(bc)
     @inline function f(i)
         copy(unpack(bc, i))
@@ -366,9 +352,9 @@ end
 end
 
 @inline function Base.copyto!(dest::ArrayPartition,
-                              bc::Broadcast.Broadcasted{ArrayPartitionStyle{Style}}) where {
-                                                                                            Style
-                                                                                            }
+    bc::Broadcast.Broadcasted{ArrayPartitionStyle{Style}}) where {
+    Style,
+}
     N = npartitions(dest, bc)
     @inline function f(i)
         copyto!(dest.x[i], unpack(bc, i))
@@ -401,15 +387,15 @@ _npartitions(args::Tuple{}) = 0
     Broadcast.Broadcasted(bc.f, unpack_args(i, bc.args))
 end
 @inline function unpack(bc::Broadcast.Broadcasted{ArrayPartitionStyle{Style}},
-                        i) where {Style}
+    i) where {Style}
     Broadcast.Broadcasted(bc.f, unpack_args(i, bc.args))
 end
 @inline function unpack(bc::Broadcast.Broadcasted{Style},
-                        i) where {Style <: Broadcast.DefaultArrayStyle}
+    i) where {Style <: Broadcast.DefaultArrayStyle}
     Broadcast.Broadcasted{Style}(bc.f, unpack_args(i, bc.args))
 end
 @inline function unpack(bc::Broadcast.Broadcasted{ArrayPartitionStyle{Style}},
-                        i) where {Style <: Broadcast.DefaultArrayStyle}
+    i) where {Style <: Broadcast.DefaultArrayStyle}
     Broadcast.Broadcasted{Style}(bc.f, unpack_args(i, bc.args))
 end
 unpack(x, ::Any) = x
@@ -438,11 +424,17 @@ function LinearAlgebra.ldiv!(A::Factorization, b::ArrayPartition)
 end
 
 @static if VERSION >= v"1.9"
-    function LinearAlgebra.ldiv!(A::LinearAlgebra.SVD{T, Tr, M}, b::ArrayPartition) where {Tr, T, M<:AbstractArray{T}}
+    function LinearAlgebra.ldiv!(A::LinearAlgebra.SVD{T, Tr, M},
+        b::ArrayPartition) where {Tr, T, M <: AbstractArray{T}}
         (x = ldiv!(A, Array(b)); copyto!(b, x))
     end
 
-    function LinearAlgebra.ldiv!(A::LinearAlgebra.QRCompactWY{T, M, C}, b::ArrayPartition) where {T<:Union{Float32, Float64, ComplexF64, ComplexF32}, M<:AbstractMatrix{T}, C<:AbstractMatrix{T}}
+    function LinearAlgebra.ldiv!(A::LinearAlgebra.QRCompactWY{T, M, C},
+        b::ArrayPartition) where {
+        T <: Union{Float32, Float64, ComplexF64, ComplexF32},
+        M <: AbstractMatrix{T},
+        C <: AbstractMatrix{T},
+    }
         (x = ldiv!(A, Array(b)); copyto!(b, x))
     end
 end
@@ -534,7 +526,7 @@ function _ldiv!(A::LowerTriangular, bb::ArrayPartition)
 end
 
 function LinearAlgebra.ldiv!(A::LowerTriangular{T, <:LinearAlgebra.Adjoint{T}},
-                             bb::ArrayPartition) where {T}
+    bb::ArrayPartition) where {T}
     _ldiv!(A, bb)
 end
 LinearAlgebra.ldiv!(A::LowerTriangular, bb::ArrayPartition) = _ldiv!(A, bb)
@@ -580,17 +572,19 @@ function LinearAlgebra.mul!(C::ArrayPartition, A::ArrayPartition, B::ArrayPartit
 end
 
 function Base.convert(::Type{ArrayPartition{T, S}},
-                      A::ArrayPartition{<:Any, <:NTuple{N, Any}}) where {N, T,
-                                                                         S <:
-                                                                         NTuple{N, Any}}
+    A::ArrayPartition{<:Any, <:NTuple{N, Any}}) where {N, T,
+    S <:
+    NTuple{N, Any}}
     return ArrayPartition{T, S}(ntuple((@inline i -> convert(S.parameters[i], A.x[i])),
-                                       Val(N)))
+        Val(N)))
 end
 
-@generated function Base.length(::Type{<:ArrayPartition{F, T}}) where {F, N,
-                                                                       T <: NTuple{N,
-                                                                              StaticArraysCore.StaticArray
-                                                                              }}
+@generated function Base.length(::Type{
+    <:ArrayPartition{F, T},
+}) where {F, N,
+    T <: NTuple{N,
+        StaticArraysCore.StaticArray,
+    }}
     sum_expr = Expr(:call, :+)
     for param in T.parameters
         push!(sum_expr.args, :(length($param)))
diff --git a/src/tabletraits.jl b/src/tabletraits.jl
index 71ba06aa..cd6e9dd8 100644
--- a/src/tabletraits.jl
+++ b/src/tabletraits.jl
@@ -7,16 +7,15 @@ function Tables.rows(A::AbstractDiffEqArray)
         N = length(A.u[1])
         names = [
             :timestamp,
-            (!(A.sc isa SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}) ?
-             (states(A.sc)[i] for i in 1:N) :
-             (Symbol("value", i) for i in 1:N))...,
+            (isempty(variable_symbols(A)) ?
+             (Symbol("value", i) for i in 1:N) :
+             Symbol.(variable_symbols(A)))...,
         ]
         types = Type[eltype(A.t), (eltype(A.u[1]) for _ in 1:N)...]
     else
         names = [
             :timestamp,
-            !(A.sc isa SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}) ?
-            states(A.sc)[1] : :value,
+            (isempty(variable_symbols(A)) ? :value : Symbol(variable_symbols(A)[1])),
         ]
         types = Type[eltype(A.t), VT]
     end
@@ -37,7 +36,7 @@ struct AbstractDiffEqArrayRows{T, U}
 end
 function AbstractDiffEqArrayRows(names, types, t, u)
     AbstractDiffEqArrayRows(Symbol.(names), types,
-                            Dict(Symbol(nm) => i for (i, nm) in enumerate(names)), t, u)
+        Dict(Symbol(nm) => i for (i, nm) in enumerate(names)), t, u)
 end
 
 Base.length(x::AbstractDiffEqArrayRows) = length(x.u)
@@ -45,7 +44,7 @@ function Base.eltype(::Type{AbstractDiffEqArrayRows{T, U}}) where {T, U}
     AbstractDiffEqArrayRow{eltype(T), eltype(U)}
 end
 function Base.iterate(x::AbstractDiffEqArrayRows,
-                      (t_state, u_state) = (iterate(x.t), iterate(x.u)))
+    (t_state, u_state) = (iterate(x.t), iterate(x.u)))
     t_state === nothing && return nothing
     u_state === nothing && return nothing
     t, _t_state = t_state
diff --git a/src/utils.jl b/src/utils.jl
index bc8d7951..18dad100 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -10,7 +10,7 @@ function recursivecopy(a)
     deepcopy(a)
 end
 function recursivecopy(a::Union{StaticArraysCore.SVector, StaticArraysCore.SMatrix,
-                                StaticArraysCore.SArray, Number})
+    StaticArraysCore.SArray, Number})
     copy(a)
 end
 function recursivecopy(a::AbstractArray{T, N}) where {T <: Number, N}
@@ -37,9 +37,9 @@ like `copy!` on arrays of scalars.
 function recursivecopy! end
 
 function recursivecopy!(b::AbstractArray{T, N},
-                        a::AbstractArray{T2, N}) where {T <: StaticArraysCore.StaticArray,
-                                                        T2 <: StaticArraysCore.StaticArray,
-                                                        N}
+    a::AbstractArray{T2, N}) where {T <: StaticArraysCore.StaticArray,
+    T2 <: StaticArraysCore.StaticArray,
+    N}
     @inbounds for i in eachindex(a)
         # TODO: Check for `setindex!`` and use `copy!(b[i],a[i])` or `b[i] = a[i]`, see #19
         b[i] = copy(a[i])
@@ -47,18 +47,18 @@ function recursivecopy!(b::AbstractArray{T, N},
 end
 
 function recursivecopy!(b::AbstractArray{T, N},
-                        a::AbstractArray{T2, N}) where {T <: Enum, T2 <: Enum, N}
+    a::AbstractArray{T2, N}) where {T <: Enum, T2 <: Enum, N}
     copyto!(b, a)
 end
 
 function recursivecopy!(b::AbstractArray{T, N},
-                        a::AbstractArray{T2, N}) where {T <: Number, T2 <: Number, N}
+    a::AbstractArray{T2, N}) where {T <: Number, T2 <: Number, N}
     copyto!(b, a)
 end
 
 function recursivecopy!(b::AbstractArray{T, N},
-                        a::AbstractArray{T2, N}) where {T <: AbstractArray,
-                                                        T2 <: AbstractArray, N}
+    a::AbstractArray{T2, N}) where {T <: AbstractArray,
+    T2 <: AbstractArray, N}
     if ArrayInterface.ismutable(T)
         @inbounds for i in eachindex(b, a)
             recursivecopy!(b[i], a[i])
@@ -79,32 +79,32 @@ A recursive `fill!` function.
 function recursivefill! end
 
 function recursivefill!(b::AbstractArray{T, N},
-                        a::T2) where {T <: StaticArraysCore.StaticArray,
-                                      T2 <: StaticArraysCore.StaticArray, N}
+    a::T2) where {T <: StaticArraysCore.StaticArray,
+    T2 <: StaticArraysCore.StaticArray, N}
     @inbounds for i in eachindex(b)
         b[i] = copy(a)
     end
 end
 
 function recursivefill!(bs::AbstractVectorOfArray{T, N},
-                        a::T2) where {T <: StaticArraysCore.StaticArray,
-                                      T2 <: StaticArraysCore.StaticArray, N}
+    a::T2) where {T <: StaticArraysCore.StaticArray,
+    T2 <: StaticArraysCore.StaticArray, N}
     @inbounds for b in bs, i in eachindex(b)
         b[i] = copy(a)
     end
 end
 
 function recursivefill!(b::AbstractArray{T, N},
-                        a::T2) where {T <: StaticArraysCore.SArray,
-                                      T2 <: Union{Number, Bool}, N}
+    a::T2) where {T <: StaticArraysCore.SArray,
+    T2 <: Union{Number, Bool}, N}
     @inbounds for i in eachindex(b)
         b[i] = fill(a, typeof(b[i]))
     end
 end
 
 function recursivefill!(bs::AbstractVectorOfArray{T, N},
-                        a::T2) where {T <: StaticArraysCore.SArray,
-                                      T2 <: Union{Number, Bool}, N}
+    a::T2) where {T <: StaticArraysCore.SArray,
+    T2 <: Union{Number, Bool}, N}
     @inbounds for b in bs, i in eachindex(b)
         b[i] = fill(a, typeof(b[i]))
     end
@@ -115,8 +115,8 @@ function recursivefill!(b::AbstractArray{T, N}, a::T2) where {T <: Enum, T2 <: E
 end
 
 function recursivefill!(b::AbstractArray{T, N},
-                        a::T2) where {T <: Union{Number, Bool}, T2 <: Union{Number, Bool}, N
-                                      }
+    a::T2) where {T <: Union{Number, Bool}, T2 <: Union{Number, Bool}, N
+}
     fill!(b, a)
 end
 
@@ -204,7 +204,7 @@ function copyat_or_push!(a::AbstractVector{T}, i::Int, x, perform_copy = true) w
 end
 
 function copyat_or_push!(a::AbstractVector{T}, i::Int, x,
-                         nc::Type{Val{perform_copy}}) where {T, perform_copy}
+    nc::Type{Val{perform_copy}}) where {T, perform_copy}
     copyat_or_push!(a, i, x, perform_copy)
 end
 
diff --git a/src/vector_of_array.jl b/src/vector_of_array.jl
index 64126f05..42560d6b 100644
--- a/src/vector_of_array.jl
+++ b/src/vector_of_array.jl
@@ -53,52 +53,71 @@ A[1, :]  # all time periods for f(t)
 A.t
 ```
 """
-mutable struct DiffEqArray{T, N, A, B, C, E, F} <: AbstractDiffEqArray{T, N, A}
+mutable struct DiffEqArray{T, N, A, B, F, S} <: AbstractDiffEqArray{T, N, A}
     u::A # A <: AbstractVector{<: AbstractArray{T, N - 1}}
     t::B
-    sc::C
-    observed::E
     p::F
+    sys::S
 end
 ### Abstract Interface
 struct AllObserved
 end
 
-# extended by Symbolcs
-issymbollike(::Any) = false
-issymbollike(::Union{Symbol, AllObserved}) = true
-
-function Base.Array(VA::AbstractVectorOfArray{T, N, A}) where {T, N,
-                                                               A <: AbstractVector{
-                                                                              <:AbstractVector
-                                                                              }}
+function Base.Array(VA::AbstractVectorOfArray{
+    T,
+    N,
+    A,
+}) where {T, N,
+    A <: AbstractVector{
+        <:AbstractVector,
+    }}
     reduce(hcat, VA.u)
 end
-function Base.Array(VA::AbstractVectorOfArray{T, N, A}) where {T, N,
-                                                               A <:
-                                                               AbstractVector{<:Number}}
+function Base.Array(VA::AbstractVectorOfArray{
+    T,
+    N,
+    A,
+}) where {T, N,
+    A <:
+    AbstractVector{<:Number}}
     VA.u
 end
-function Base.Matrix(VA::AbstractVectorOfArray{T, N, A}) where {T, N,
-                                                                A <: AbstractVector{
-                                                                               <:AbstractVector
-                                                                               }}
+function Base.Matrix(VA::AbstractVectorOfArray{
+    T,
+    N,
+    A,
+}) where {T, N,
+    A <: AbstractVector{
+        <:AbstractVector,
+    }}
     reduce(hcat, VA.u)
 end
-function Base.Matrix(VA::AbstractVectorOfArray{T, N, A}) where {T, N,
-                                                                A <:
-                                                                AbstractVector{<:Number}}
+function Base.Matrix(VA::AbstractVectorOfArray{
+    T,
+    N,
+    A,
+}) where {T, N,
+    A <:
+    AbstractVector{<:Number}}
     Matrix(VA.u)
 end
-function Base.Vector(VA::AbstractVectorOfArray{T, N, A}) where {T, N,
-                                                                A <: AbstractVector{
-                                                                               <:AbstractVector
-                                                                               }}
+function Base.Vector(VA::AbstractVectorOfArray{
+    T,
+    N,
+    A,
+}) where {T, N,
+    A <: AbstractVector{
+        <:AbstractVector,
+    }}
     vec(reduce(hcat, VA.u))
 end
-function Base.Vector(VA::AbstractVectorOfArray{T, N, A}) where {T, N,
-                                                                A <:
-                                                                AbstractVector{<:Number}}
+function Base.Vector(VA::AbstractVectorOfArray{
+    T,
+    N,
+    A,
+}) where {T, N,
+    A <:
+    AbstractVector{<:Number}}
     VA.u
 end
 function Base.Array(VA::AbstractVectorOfArray)
@@ -119,198 +138,255 @@ function VectorOfArray(vec::AbstractVector{VT}) where {T, N, VT <: AbstractArray
     VectorOfArray{T, N + 1, typeof(vec)}(vec)
 end
 
-function DiffEqArray(vec::AbstractVector{T}, ts, ::NTuple{N, Int}, syms = nothing,
-                     indepsym = nothing, observed = nothing, p = nothing) where {T, N}
-    sc = if isnothing(indepsym) || indepsym isa AbstractArray
-        SymbolCache{typeof(syms), typeof(indepsym), Nothing}(syms, indepsym, nothing)
-    else
-        SymbolCache{typeof(syms), Vector{typeof(indepsym)}, Nothing}(syms, [indepsym],
-                                                                     nothing)
-    end
-    DiffEqArray{eltype(T), N, typeof(vec), typeof(ts), typeof(sc), typeof(observed),
-                typeof(p)}(vec, ts, sc, observed, p)
+function DiffEqArray(vec::AbstractVector{T},
+    ts,
+    ::NTuple{N, Int},
+    p = nothing,
+    sys = nothing) where {T, N}
+    DiffEqArray{eltype(T), N, typeof(vec), typeof(ts), typeof(p), typeof(sys)}(vec,
+        ts,
+        p,
+        sys)
 end
 # Assume that the first element is representative of all other elements
-function DiffEqArray(vec::AbstractVector, ts::AbstractVector, syms = nothing,
-                     indepsym = nothing, observed = nothing, p = nothing)
-    DiffEqArray(vec, ts, (size(vec[1])..., length(vec)), syms, indepsym, observed, p)
-end
-function DiffEqArray(vec::AbstractVector{VT}, ts::AbstractVector, syms = nothing,
-                     indepsym = nothing, observed = nothing,
-                     p = nothing) where {T, N, VT <: AbstractArray{T, N}}
-    sc = if isnothing(indepsym) || indepsym isa AbstractArray
-        SymbolCache{typeof(syms), typeof(indepsym), Nothing}(syms, indepsym, nothing)
-    else
-        SymbolCache{typeof(syms), Vector{typeof(indepsym)}, Nothing}(syms, [indepsym],
-                                                                     nothing)
-    end
-    DiffEqArray{T, N + 1, typeof(vec), typeof(ts), typeof(sc), typeof(observed), typeof(p)}(vec,
-                                                                                            ts,
-                                                                                            sc,
-                                                                                            observed,
-                                                                                            p)
-end
 
-# Interface for the linear indexing. This is just a view of the underlying nested structure
-@inline Base.firstindex(VA::AbstractVectorOfArray) = firstindex(VA.u)
-@inline Base.lastindex(VA::AbstractVectorOfArray) = lastindex(VA.u)
+function DiffEqArray(vec::AbstractVector,
+    ts::AbstractVector,
+    p = nothing,
+    sys = nothing;
+    variables = nothing,
+    parameters = nothing,
+    independent_variables = nothing)
+    sys = something(sys,
+        SymbolCache(something(variables, []),
+            something(parameters, []),
+            something(independent_variables, [])))
+    _size = size(vec[1])
+    return DiffEqArray{
+        eltype(eltype(vec)),
+        length(_size),
+        typeof(vec),
+        typeof(ts),
+        typeof(p),
+        typeof(sys),
+    }(vec,
+        ts,
+        p,
+        sys)
+end
+
+function DiffEqArray(vec::AbstractVector{VT},
+    ts::AbstractVector,
+    p = nothing,
+    sys = nothing;
+    variables = nothing,
+    parameters = nothing,
+    independent_variables = nothing) where {T, N, VT <: AbstractArray{T, N}}
+    sys = something(sys, SymbolCache(something(variables, []),
+        something(parameters, []),
+        something(independent_variables, [])))
+    return DiffEqArray{
+        eltype(eltype(vec)),
+        N + 1,
+        typeof(vec),
+        typeof(ts),
+        typeof(p),
+        typeof(sys),
+    }(vec,
+        ts,
+        p,
+        sys)
+end
+
+SymbolicIndexingInterface.parameter_values(A::AbstractDiffEqArray) = A.p
+SymbolicIndexingInterface.symbolic_container(A::AbstractDiffEqArray) = A.sys
+
+Base.IndexStyle(A::AbstractVectorOfArray) = Base.IndexStyle(typeof(A))
+Base.IndexStyle(::Type{<:AbstractVectorOfArray}) = IndexCartesian()
 
 @inline Base.length(VA::AbstractVectorOfArray) = length(VA.u)
-@inline Base.eachindex(VA::AbstractVectorOfArray) = Base.OneTo(length(VA.u))
-@inline Base.IteratorSize(VA::AbstractVectorOfArray) = Base.HasLength()
-# Linear indexing will be over the container elements, not the individual elements
-# unlike an true AbstractArray
-Base.@propagate_inbounds function Base.getindex(VA::AbstractVectorOfArray{T, N},
-                                                I::Int) where {T, N}
-    VA.u[I]
+@inline function Base.eachindex(VA::AbstractVectorOfArray)
+    return eachindex(VA.u)
 end
-Base.@propagate_inbounds function Base.getindex(VA::AbstractVectorOfArray{T, N},
-                                                I::Colon) where {T, N}
-    VA.u[I]
+@inline Base.IteratorSize(::Type{<:AbstractVectorOfArray}) = Base.HasLength()
+@inline Base.first(VA::AbstractVectorOfArray) = first(VA.u)
+@inline Base.last(VA::AbstractVectorOfArray) = last(VA.u)
+function Base.firstindex(VA::AbstractVectorOfArray)
+    Base.depwarn(Linear indexing of `AbstractVectorOfArray` is deprecated. Change `A[i]` to `A.u[i]` ",, :firstindex)
+    return firstindex(VA.u)
 end
-Base.@propagate_inbounds function Base.getindex(VA::AbstractDiffEqArray{T, N},
-                                                I::Colon) where {T, N}
-    VA.u[I]
-end
-Base.@propagate_inbounds function Base.getindex(VA::AbstractVectorOfArray{T, N},
-                                                I::AbstractArray{Int}) where {T, N}
-    VectorOfArray(VA.u[I])
-end
-Base.@propagate_inbounds function Base.getindex(VA::AbstractDiffEqArray{T, N},
-                                                I::AbstractArray{Int}) where {T, N}
-    DiffEqArray(VA.u[I], VA.t[I])
-end
-Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray{T, N},
-                                                I::Union{Int, AbstractArray{Int},
-                                                         CartesianIndex, Colon, BitArray,
-                                                         AbstractArray{Bool}}...) where {T,
-                                                                                         N}
-    RecursiveArrayTools.VectorOfArray(A.u)[I...]
+
+function Base.lastindex(VA::AbstractVectorOfArray)
+    Base.depwarn("Linear indexing of `AbstractVectorOfArray` is deprecated. Change `A[i]` to `A.u[i]` ", :lastindex)
+    return lastindex(VA.u)
 end
 
+@deprecate Base.getindex(A::AbstractVectorOfArray, I::Int) Base.getindex(A, :, I) false
+
+@deprecate Base.getindex(A::AbstractVectorOfArray, I::AbstractArray{Int}) Base.getindex(A, :, I) false
+
+@deprecate Base.getindex(A::AbstractDiffEqArray, I::AbstractArray{Int}) Base.getindex(A, :, I) false
+
+@deprecate Base.getindex(A::AbstractDiffEqArray, i::Int) Base.getindex(A, :, i) false
+
 __parameterless_type(T) = Base.typename(T).wrapper
 Base.@propagate_inbounds function Base.getindex(A::AbstractVectorOfArray{T, N},
-                                                I::Colon...) where {T, N}
+    ::NotSymbolic, I::Colon...) where {T, N}
     @assert length(I) == ndims(A.u[1]) + 1
     vecs = vec.(A.u)
     return Adapt.adapt(__parameterless_type(T),
-                       reshape(reduce(hcat, vecs), size(A.u[1])..., length(A.u)))
+        reshape(reduce(hcat, vecs), size(A.u[1])..., length(A.u)))
 end
 
 Base.@propagate_inbounds function Base.getindex(A::AbstractVectorOfArray{T, N},
-                                                I::AbstractArray{Bool},
-                                                J::Colon...) where {T, N}
+    ::NotSymbolic, I::AbstractArray{Bool},
+    J::Colon...) where {T, N}
     @assert length(J) == ndims(A.u[1]) + 1 - ndims(I)
     @assert size(I) == size(A)[1:(ndims(A) - length(J))]
     return A[ntuple(x -> Colon(), ndims(A))...][I, J...]
 end
 
-Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray{T, N}, i::Int,
-                                                ::Colon) where {T, N}
-    [A.u[j][i] for j in 1:length(A)]
-end
-Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray{T, N}, ::Colon,
-                                                i::Int) where {T, N}
-    A.u[i]
-end
-Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray{T, N}, i::Int,
-                                                II::AbstractArray{Int}) where {T, N}
-    [A.u[j][i] for j in II]
-end
-Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray{T, N},
-                                                sym) where {T, N}
-    if issymbollike(sym) && !isnothing(A.sc)
-        if is_indep_sym(A.sc, sym)
-            return A.t
-        elseif is_state_sym(A.sc, sym)
-            return getindex.(A.u, state_sym_to_index(A.sc, sym))
-        elseif is_param_sym(A.sc, sym)
-            return A.p[param_sym_to_index(A.sc, sym)]
-        elseif A.observed !== nothing
-            return observed(A, sym, :)
-        end
-    elseif all(issymbollike, sym) && !isnothing(A.sc)
-        if all(Base.Fix1(is_param_sym, A.sc), sym)
-            return getindex.((A,), sym)
-        else
-            return [getindex.((A,), sym, i) for i in eachindex(A.t)]
-        end
-    end
-    return getindex.(A.u, sym)
-end
-Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray{T, N}, sym,
-                                                args...) where {T, N}
-    if issymbollike(sym) && !isnothing(A.sc)
-        if is_indep_sym(A.sc, sym)
-            return A.t[args...]
-        elseif is_state_sym(A.sc, sym)
-            return A[sym][args...]
-        elseif A.observed !== nothing
-            return observed(A, sym, args...)
-        end
-    elseif all(issymbollike, sym) && !isnothing(A.sc)
-        return reduce(vcat, map(s -> A[s, args...]', sym))
-    end
-    return getindex.(A.u, sym)
-end
-Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray{T, N},
-                                                I::Int...) where {T, N}
-    A.u[I[end]][Base.front(I)...]
+# Need two of each methods to avoid ambiguities
+Base.@propagate_inbounds function Base.getindex(A::AbstractVectorOfArray, ::NotSymbolic, ::Colon, I::Int)
+    A.u[I]
 end
-Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray{T, N},
-                                                i::Int) where {T, N}
-    A.u[i]
+
+Base.@propagate_inbounds function Base.getindex(A::AbstractVectorOfArray, ::NotSymbolic, I::Union{Int,AbstractArray{Int},AbstractArray{Bool},Colon}...)
+    if last(I) isa Int
+        A.u[last(I)][Base.front(I)...]
+    else
+        stack(getindex.(A.u[last(I)], tuple.(Base.front(I))...))
+    end
 end
-Base.@propagate_inbounds function Base.getindex(VA::AbstractDiffEqArray{T, N},
-                                                ii::CartesianIndex) where {T, N}
+Base.@propagate_inbounds function Base.getindex(VA::AbstractVectorOfArray, ::NotSymbolic, ii::CartesianIndex)
     ti = Tuple(ii)
     i = last(ti)
     jj = CartesianIndex(Base.front(ti))
     return VA.u[i][jj]
 end
 
-function observed(A::AbstractDiffEqArray{T, N}, sym, i::Int) where {T, N}
-    A.observed(sym, A.u[i], A.p, A.t[i])
+Base.@propagate_inbounds function Base.getindex(A::AbstractVectorOfArray, ::NotSymbolic, ::Colon, I::Union{AbstractArray{Int},AbstractArray{Bool}})
+    VectorOfArray(A.u[I])
 end
-function observed(A::AbstractDiffEqArray{T, N}, sym, i::AbstractArray{Int}) where {T, N}
-    A.observed.((sym,), A.u[i], (A.p,), A.t[i])
+
+Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray, ::NotSymbolic, ::Colon, I::Union{AbstractArray{Int},AbstractArray{Bool}})
+    DiffEqArray(A.u[I], A.t[I], parameter_values(A), symbolic_container(A))
 end
-function observed(A::AbstractDiffEqArray{T, N}, sym, ::Colon) where {T, N}
-    A.observed.((sym,), A.u, (A.p,), A.t)
+
+# Symbolic Indexing Methods
+Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray, ::ScalarSymbolic, sym)
+    if is_independent_variable(A, sym)
+        return A.t
+    elseif is_variable(A, sym)
+        if constant_structure(A)
+            return getindex.(A.u, variable_index(A, sym))
+        else
+            return getindex.(A.u, variable_index.((A,), (sym,), eachindex(A.t)))
+        end
+    elseif is_parameter(A, sym)
+        Base.depwarn("Indexing with parameters is deprecated. Use `getp(A, $sym)` for parameter indexing.", :parameter_getindex)
+        return getp(A, sym)(A)
+    elseif is_observed(A, sym)
+        return observed(A, sym).(A.u, (parameter_values(A),), A.t)
+    else
+        # NOTE: this is basically just for LabelledArrays. It's better if this
+        # were an error. Should we make an extension for LabelledArrays handling
+        # this case?
+        return getindex.(A.u, sym)
+    end
 end
 
-Base.@propagate_inbounds function Base.getindex(VA::AbstractVectorOfArray{T, N}, i::Int,
-                                                ::Colon) where {T, N}
-    [VA.u[j][i] for j in 1:length(VA)]
+Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray, ::ScalarSymbolic, sym, args...)
+    if is_independent_variable(A, sym)
+        return A.t[args...]
+    elseif is_variable(A, sym)
+        return A[sym][args...]
+    elseif is_observed(A, sym)
+        u = A.u[args...]
+        t = A.t[args...]
+        observed_fn = observed(A, sym)
+        if t isa AbstractArray
+            return observed_fn.(u, (parameter_values(A),), t)
+        else
+            return observed_fn(u, parameter_values(A), t)
+        end
+    else
+        # NOTE: this is basically just for LabelledArrays. It's better if this
+        # were an error. Should we make an extension for LabelledArrays handling
+        # this case?
+        return getindex.(A.u[args...], sym)
+    end
 end
-Base.@propagate_inbounds function Base.getindex(VA::AbstractVectorOfArray{T, N},
-                                                ii::CartesianIndex) where {T, N}
-    ti = Tuple(ii)
-    i = last(ti)
-    jj = CartesianIndex(Base.front(ti))
-    return VA.u[i][jj]
+
+
+Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray, ::ArraySymbolic, sym, args...)
+    return getindex(A, collect(sym), args...)
+end
+
+Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray, ::ScalarSymbolic, sym::Union{Tuple,AbstractArray})
+    if all(x -> is_parameter(A, x), sym)
+        Base.depwarn("Indexing with parameters is deprecated. Use `getp(A, $sym)` for parameter indexing.", :parameter_getindex)
+        return getp(A, sym)(A)
+    else
+        return [getindex.((A,), sym, i) for i in eachindex(A.t)]
+    end
+end
+
+Base.@propagate_inbounds function Base.getindex(A::AbstractDiffEqArray, ::ScalarSymbolic, sym::Union{Tuple,AbstractArray}, args...)
+    return reduce(vcat, map(s -> A[s, args...]', sym))
+end
+
+Base.@propagate_inbounds function Base.getindex(A::AbstractVectorOfArray, _arg, args...)
+    symtype = symbolic_type(_arg)
+    elsymtype = symbolic_type(eltype(_arg))
+
+    if symtype != NotSymbolic()
+        return Base.getindex(A, symtype, _arg, args...)
+    else
+        return Base.getindex(A, elsymtype, _arg, args...)
+    end
+end
+
+function _observed(A::AbstractDiffEqArray{T, N}, sym, i::Int) where {T, N}
+    observed(A, sym)(A.u[i], A.p, A.t[i])
 end
+function _observed(A::AbstractDiffEqArray{T, N}, sym, i::AbstractArray{Int}) where {T, N}
+    observed(A, sym).(A.u[i], (A.p,), A.t[i])
+end
+function _observed(A::AbstractDiffEqArray{T, N}, sym, ::Colon) where {T, N}
+    observed(A, sym).(A.u, (A.p,), A.t)
+end
+
 Base.@propagate_inbounds function Base.setindex!(VA::AbstractVectorOfArray{T, N}, v,
-                                                 I::Int) where {T, N}
+    ::Colon, I::Int) where {T, N}
     VA.u[I] = v
 end
+
+@deprecate Base.setindex!(VA::AbstractVectorOfArray, v, I::Int) Base.setindex!(VA, v, :, I) false
+
 Base.@propagate_inbounds function Base.setindex!(VA::AbstractVectorOfArray{T, N}, v,
-                                                 I::Colon) where {T, N}
+    ::Colon, I::Colon) where {T, N}
     VA.u[I] = v
 end
+
+@deprecate Base.setindex!(VA::AbstractVectorOfArray, v, I::Colon) Base.setindex!(VA, v, :, I) false
+
 Base.@propagate_inbounds function Base.setindex!(VA::AbstractVectorOfArray{T, N}, v,
-                                                 I::AbstractArray{Int}) where {T, N}
+    ::Colon, I::AbstractArray{Int}) where {T, N}
     VA.u[I] = v
 end
+
+@deprecate Base.setindex!(VA::AbstractVectorOfArray, v, I::AbstractArray{Int}) Base.setindex!(VA, v, :, I) false
+
 Base.@propagate_inbounds function Base.setindex!(VA::AbstractVectorOfArray{T, N}, v, i::Int,
-                                                 ::Colon) where {T, N}
-    for j in 1:length(VA)
+    ::Colon) where {T, N}
+    for j in 1:length(VA.u)
         VA.u[j][i] = v[j]
     end
     return v
 end
 Base.@propagate_inbounds function Base.setindex!(VA::AbstractVectorOfArray{T, N}, x,
-                                                 ii::CartesianIndex) where {T, N}
+    ii::CartesianIndex) where {T, N}
     ti = Tuple(ii)
     i = last(ti)
     jj = CartesianIndex(Base.front(ti))
@@ -319,91 +395,164 @@ end
 
 # Interface for the two-dimensional indexing, a more standard AbstractArray interface
 @inline Base.size(VA::AbstractVectorOfArray) = (size(VA.u[1])..., length(VA.u))
-Base.@propagate_inbounds function Base.getindex(VA::AbstractVectorOfArray{T, N},
-                                                I::Int...) where {T, N}
-    VA.u[I[end]][Base.front(I)...]
-end
-Base.@propagate_inbounds function Base.getindex(VA::AbstractVectorOfArray{T, N}, ::Colon,
-                                                I::Int) where {T, N}
-    VA.u[I]
-end
+Base.axes(VA::AbstractVectorOfArray) = Base.OneTo.(size(VA))
+Base.axes(VA::AbstractVectorOfArray, d::Int) = Base.OneTo(size(VA)[d])
+
 Base.@propagate_inbounds function Base.setindex!(VA::AbstractVectorOfArray{T, N}, v,
-                                                 I::Int...) where {T, N}
+    I::Int...) where {T, N}
     VA.u[I[end]][Base.front(I)...] = v
 end
 
+function Base.:(==)(A::AbstractVectorOfArray, B::AbstractVectorOfArray)
+    return A.u == B.u
+end
+function Base.:(==)(A::AbstractVectorOfArray, B::AbstractArray)
+    return A.u == B
+end
+Base.:(==)(A::AbstractArray, B::AbstractVectorOfArray) = B == A
+
 # The iterator will be over the subarrays of the container, not the individual elements
 # unlike an true AbstractArray
 function Base.iterate(VA::AbstractVectorOfArray, state = 1)
-    state >= length(VA.u) + 1 ? nothing : (VA[state], state + 1)
+    state >= length(VA.u) + 1 ? nothing : (VA[:, state], state + 1)
 end
 tuples(VA::DiffEqArray) = tuple.(VA.t, VA.u)
 
 # Growing the array simply adds to the container vector
 function Base.copy(VA::AbstractDiffEqArray)
     typeof(VA)(copy(VA.u),
-               copy(VA.t),
-               (VA.sc === nothing) ? nothing : copy(VA.sc),
-               (VA.observed === nothing) ? nothing : copy(VA.observed),
-               (VA.p === nothing) ? nothing : copy(VA.p))
+        copy(VA.t),
+        (VA.p === nothing) ? nothing : copy(VA.p),
+        (VA.sys === nothing) ? nothing : copy(VA.sys))
 end
 Base.copy(VA::AbstractVectorOfArray) = typeof(VA)(copy(VA.u))
 Base.sizehint!(VA::AbstractVectorOfArray{T, N}, i) where {T, N} = sizehint!(VA.u, i)
 
 Base.reverse!(VA::AbstractVectorOfArray) = reverse!(VA.u)
 Base.reverse(VA::VectorOfArray) = VectorOfArray(reverse(VA.u))
-Base.reverse(VA::DiffEqArray) = DiffEqArray(reverse(VA.u), VA.t, VA.sc, VA.observed, VA.p)
+Base.reverse(VA::DiffEqArray) = DiffEqArray(reverse(VA.u), VA.t, VA.p, VA.sys)
 
 function Base.push!(VA::AbstractVectorOfArray{T, N}, new_item::AbstractArray) where {T, N}
     push!(VA.u, new_item)
 end
 
 function Base.append!(VA::AbstractVectorOfArray{T, N},
-                      new_item::AbstractVectorOfArray{T, N}) where {T, N}
+    new_item::AbstractVectorOfArray{T, N}) where {T, N}
     for item in copy(new_item)
         push!(VA, item)
     end
     return VA
 end
 
+# AbstractArray methods
+function Base.view(A::AbstractVectorOfArray, I::Vararg{Any,M}) where {M}
+    @inline
+    J = map(i->Base.unalias(A,i), to_indices(A, I))
+    @boundscheck checkbounds(A, J...)
+    SubArray(IndexStyle(A), A, J, Base.index_dimsum(J...))
+end
+Base.check_parent_index_match(::RecursiveArrayTools.AbstractVectorOfArray{T,N}, ::NTuple{N,Bool}) where {T,N} = nothing
+Base.ndims(::AbstractVectorOfArray{T, N}) where {T, N} = N
+function Base.checkbounds(::Type{Bool}, VA::AbstractVectorOfArray, idx...)
+    if checkbounds(Bool, VA.u, last(idx))
+        if last(idx) isa Integer
+            return all(checkbounds.(Bool, (VA.u[last(idx)],), Base.front(idx)))
+        else
+            return all(checkbounds.(Bool, VA.u[last(idx)], Base.front(idx)))
+        end
+    end
+    return false
+end
+function Base.checkbounds(VA::AbstractVectorOfArray, idx...)
+    checkbounds(Bool, VA, idx...) || throw(BoundsError(VA, idx))
+end
+
+# Operations
+function Base.isapprox(A::AbstractVectorOfArray,
+    B::Union{AbstractVectorOfArray, AbstractArray};
+    kwargs...)
+    return all(isapprox.(A, B; kwargs...))
+end
+
+function Base.isapprox(A::AbstractArray, B::AbstractVectorOfArray; kwargs...)
+    return all(isapprox.(A, B; kwargs...))
+end
+
+for op in [:(Base.:-), :(Base.:+)]
+    @eval function ($op)(A::AbstractVectorOfArray, B::AbstractVectorOfArray)
+        ($op).(A, B)
+    end
+    @eval Base.@propagate_inbounds function ($op)(A::AbstractVectorOfArray,
+        B::AbstractArray)
+        @boundscheck length(A) == length(B)
+        VectorOfArray([($op).(a, b) for (a, b) in zip(A, B)])
+    end
+    @eval Base.@propagate_inbounds function ($op)(A::AbstractArray, B::AbstractVectorOfArray)
+        @boundscheck length(A) == length(B)
+        VectorOfArray([($op).(a, b) for (a, b) in zip(A, B)])
+    end
+end
+
+for op in [:(Base.:/), :(Base.:\), :(Base.:*)]
+    if op !== :(Base.:/)
+        @eval ($op)(A::Number, B::AbstractVectorOfArray) = ($op).(A, B)
+    end
+    if op !== :(Base.:\)
+        @eval ($op)(A::AbstractVectorOfArray, B::Number) = ($op).(A, B)
+    end
+end
+
+function Base.CartesianIndices(VA::AbstractVectorOfArray)
+    if !allequal(size.(VA.u))
+        error("CartesianIndices only valid for non-ragged arrays")
+    end
+    return CartesianIndices((size(VA.u[1])..., length(VA.u)))
+end
+
 # Tools for creating similar objects
+Base.eltype(::VectorOfArray{T}) where {T} = T
+# TODO: Is there a better way to do this?
+@inline function Base.similar(VA::AbstractVectorOfArray, args...)
+    if args[end] isa Type
+        return Base.similar(eltype(VA)[], args..., size(VA))
+    else
+        return Base.similar(eltype(VA)[], args...)
+    end
+end
 @inline function Base.similar(VA::VectorOfArray, ::Type{T} = eltype(VA)) where {T}
-    VectorOfArray([similar(VA[i], T) for i in eachindex(VA)])
+    VectorOfArray([similar(VA[:, i], T) for i in eachindex(VA.u)])
 end
 recursivecopy(VA::VectorOfArray) = VectorOfArray(copy.(VA.u))
 
 # fill!
 # For DiffEqArray it ignores ts and fills only u
 function Base.fill!(VA::AbstractVectorOfArray, x)
-    for i in eachindex(VA)
-        if VA[i] isa AbstractArray
-            fill!(VA[i], x)
+    for i in 1:length(VA.u)
+        if VA[:, i] isa AbstractArray
+            fill!(VA[:, i], x)
         else
-            VA[i] = x
+            VA[:, i] = x
         end
     end
     return VA
 end
 
-function Base._reshape(parent::VectorOfArray, dims::Base.Dims)
-    n = prod(size(parent))
-    prod(dims) == n || Base._throw_dmrs(n, "size", dims)
-    Base.__reshape((parent, IndexStyle(parent)), dims)
-end
+Base.reshape(A::VectorOfArray, dims...) = Base.reshape(Array(A), dims...)
 
 # Need this for ODE_DEFAULT_UNSTABLE_CHECK from DiffEqBase to work properly
-@inline Base.any(f, VA::AbstractVectorOfArray) = any(any(f, VA[i]) for i in eachindex(VA))
-@inline Base.all(f, VA::AbstractVectorOfArray) = all(all(f, VA[i]) for i in eachindex(VA))
-@inline function Base.any(f::Function, VA::AbstractVectorOfArray)
-    any(any(f, VA[i]) for i in eachindex(VA))
-end
-@inline function Base.all(f::Function, VA::AbstractVectorOfArray)
-    all(all(f, VA[i]) for i in eachindex(VA))
-end
+@inline Base.any(f, VA::AbstractVectorOfArray) = any(any(f, u) for u in VA.u)
+@inline Base.all(f, VA::AbstractVectorOfArray) = all(all(f, u) for u in VA.u)
 
 # conversion tools
 vecarr_to_vectors(VA::AbstractVectorOfArray) = [VA[i, :] for i in eachindex(VA[1])]
 Base.vec(VA::AbstractVectorOfArray) = vec(convert(Array, VA)) # Allocates
+# stack non-ragged arrays to convert them
+function Base.convert(::Type{Array}, VA::AbstractVectorOfArray)
+    if !allequal(size.(VA.u))
+        error("Can only convert non-ragged VectorOfArray to Array")
+    end
+    return stack(VA.u)
+end
 
 # statistics
 @inline Base.sum(f, VA::AbstractVectorOfArray) = sum(f, Array(VA))
@@ -424,8 +573,8 @@ end
 function Base.show(io::IO, m::MIME"text/plain", x::AbstractVectorOfArray)
     (println(io, summary(x), ':'); show(io, m, x.u))
 end
-function Base.summary(A::AbstractVectorOfArray)
-    string("VectorOfArray{", eltype(A), ",", ndims(A), "}")
+function Base.summary(A::AbstractVectorOfArray{T, N}) where {T, N}
+    string("VectorOfArray{", T, ",", N, "}")
 end
 
 function Base.show(io::IO, m::MIME"text/plain", x::AbstractDiffEqArray)
@@ -437,8 +586,10 @@ end
     convert(Array, VA)
 end
 @recipe function f(VA::AbstractDiffEqArray)
-    xguide --> ((VA.sc.indepsym !== nothing) ? string(VA.sc.indepsym) : "")
-    label --> ((VA.sc.syms !== nothing) ? reshape(string.(VA.sc.syms), 1, :) : "")
+    xguide --> isempty(independent_variable_symbols(VA)) ? "" :
+    independent_variable_symbols(VA)[1]
+    label --> isempty(variable_symbols(VA)) ? "" :
+    reshape(string.(variable_symbols(VA)), 1, :)
     VA.t, VA'
 end
 @recipe function f(VA::DiffEqArray{T, 1}) where {T}
@@ -458,38 +609,40 @@ VectorOfArrayStyle(::Val{N}) where {N} = VectorOfArrayStyle{N}()
 # The order is important here. We want to override Base.Broadcast.DefaultArrayStyle to return another Base.Broadcast.DefaultArrayStyle.
 Broadcast.BroadcastStyle(a::VectorOfArrayStyle, ::Base.Broadcast.DefaultArrayStyle{0}) = a
 function Broadcast.BroadcastStyle(::VectorOfArrayStyle{N},
-                                  a::Base.Broadcast.DefaultArrayStyle{M}) where {M, N}
+    a::Base.Broadcast.DefaultArrayStyle{M}) where {M, N}
     Base.Broadcast.DefaultArrayStyle(Val(max(M, N)))
 end
 function Broadcast.BroadcastStyle(::VectorOfArrayStyle{N},
-                                  a::Base.Broadcast.AbstractArrayStyle{M}) where {M, N}
+    a::Base.Broadcast.AbstractArrayStyle{M}) where {M, N}
     typeof(a)(Val(max(M, N)))
 end
 function Broadcast.BroadcastStyle(::VectorOfArrayStyle{M},
-                                  ::VectorOfArrayStyle{N}) where {M, N}
+    ::VectorOfArrayStyle{N}) where {M, N}
     VectorOfArrayStyle(Val(max(M, N)))
 end
 function Broadcast.BroadcastStyle(::Type{<:AbstractVectorOfArray{T, N}}) where {T, N}
     VectorOfArrayStyle{N}()
 end
+# make vectorofarrays broadcastable so they aren't collected
+Broadcast.broadcastable(x::AbstractVectorOfArray) = x
 
 @inline function Base.copy(bc::Broadcast.Broadcasted{<:VectorOfArrayStyle})
     bc = Broadcast.flatten(bc)
     N = narrays(bc)
     VectorOfArray(map(1:N) do i
-                      copy(unpack_voa(bc, i))
-                  end)
+        copy(unpack_voa(bc, i))
+    end)
 end
 
 @inline function Base.copyto!(dest::AbstractVectorOfArray,
-                              bc::Broadcast.Broadcasted{<:VectorOfArrayStyle})
+    bc::Broadcast.Broadcasted{<:VectorOfArrayStyle})
     bc = Broadcast.flatten(bc)
     N = narrays(bc)
     @inbounds for i in 1:N
-        if dest[i] isa AbstractArray
-            copyto!(dest[i], unpack_voa(bc, i))
+        if dest[:, i] isa AbstractArray
+            copyto!(dest[:, i], unpack_voa(bc, i))
         else
-            dest[i] = copy(unpack_voa(bc, i))
+            dest[:, i] = copy(unpack_voa(bc, i))
         end
     end
     dest
@@ -503,7 +656,7 @@ end
 Retrieve number of arrays in the AbstractVectorOfArrays of a broadcast.
 """
 narrays(A) = 0
-narrays(A::AbstractVectorOfArray) = length(A)
+narrays(A::AbstractVectorOfArray) = length(A.u)
 narrays(bc::Broadcast.Broadcasted) = _narrays(bc.args)
 narrays(A, Bs...) = common_length(narrays(A), _narrays(Bs))
 
@@ -514,7 +667,7 @@ function common_length(a, b)
       throw(DimensionMismatch("number of arrays must be equal"))))
 end
 
-_narrays(args::AbstractVectorOfArray) = length(args)
+_narrays(args::AbstractVectorOfArray) = length(args.u)
 @inline _narrays(args::Tuple) = common_length(narrays(args[1]), _narrays(Base.tail(args)))
 _narrays(args::Tuple{Any}) = _narrays(args[1])
 _narrays(::Any) = 0
diff --git a/test/basic_indexing.jl b/test/basic_indexing.jl
index 5ebd7e1e..7a5acbe6 100644
--- a/test/basic_indexing.jl
+++ b/test/basic_indexing.jl
@@ -16,8 +16,8 @@ mulX .= sqrt.(abs.(testva .* X))
 @test mulX == ref
 
 @test Array(testva) == [1 4 7
-                        2 5 8
-                        3 6 9]
+    2 5 8
+    3 6 9]
 
 @test testa[1:2, 1:2] == [1 4; 2 5]
 @test testva[1:2, 1:2] == [1 4; 2 5]
@@ -26,8 +26,8 @@ mulX .= sqrt.(abs.(testva .* X))
 t = [1, 2, 3]
 diffeq = DiffEqArray(recs, t)
 @test Array(diffeq) == [1 4 7
-                        2 5 8
-                        3 6 9]
+    2 5 8
+    3 6 9]
 @test diffeq[1:2, 1:2] == [1 4; 2 5]
 
 # # ndims == 2
@@ -36,17 +36,44 @@ recs = [rand(8) for i in 1:10]
 testa = cat(recs..., dims = 2)
 testva = VectorOfArray(recs)
 
+# Array functions
+@test size(testva) == (8, 10)
+@test axes(testva) == Base.OneTo.((8, 10))
+@test ndims(testva) == 2
+@test eltype(testva) == eltype(eltype(recs))
+testvasim = similar(testva)
+@test size(testvasim) == size(testva)
+@test eltype(testvasim) == eltype(testva)
+testvasim = similar(testva, Float32)
+@test size(testvasim) == size(testva)
+@test eltype(testvasim) == Float32
+testvb = deepcopy(testva)
+@test testva == testvb == recs
+
+# Math operations
+@test testva + testvb == testva + recs == 2testva == 2 .* recs
+@test testva - testvb == testva - recs == 0 .* recs
+@test testva / 2 == recs ./ 2
+@test 2 .\ testva == 2 .\ recs
+
 # ## Linear indexing
-@test testva[1] == testa[:, 1]
-@test testva[:] == recs
-@test testva[end] == testa[:, end]
-@test testva[2:end] == VectorOfArray([recs[i] for i in 2:length(recs)])
+@test_deprecated testva[1]
+@test_deprecated testva[1:2]
+@test_deprecated testva[begin]
+@test_deprecated testva[end]
+@test testva[begin] == testva[:, begin] == first(testva)
+@test testva[end] == testva[:, end] == last(testva)
+@test testa[:, 1] == recs[1]
+@test testva.u == recs
+@test testva[: ,2:end] == VectorOfArray([recs[i] for i in 2:length(recs)])
 
 diffeq = DiffEqArray(recs, t)
-@test diffeq[1] == testa[:, 1]
-@test diffeq[:] == recs
-@test diffeq[end] == testa[:, end]
-@test diffeq[2:end] == DiffEqArray([recs[i] for i in 2:length(recs)], t)
+@test_deprecated diffeq[1]
+@test_deprecated diffeq[1:2]
+@test diffeq[:, 1] == testa[:, 1]
+@test diffeq.u == recs
+@test diffeq[:, end] == testa[:, end]
+@test diffeq[:, 2:end] == DiffEqArray([recs[i] for i in 2:length(recs)], t)
 
 # ## (Int, Int)
 @test testa[5, 4] == testva[5, 4]
@@ -139,10 +166,10 @@ v[CartesianIndex((2, 3, 2, 3))] = 1
 v = VectorOfArray([rand(20), rand(10, 10), rand(3, 3, 3)])
 w = v .* v
 @test w isa VectorOfArray
-@test w[1] isa Vector
-@test w[1] == v[1] .* v[1]
-@test w[2] == v[2] .* v[2]
-@test w[3] == v[3] .* v[3]
+@test w[:, 1] isa Vector
+@test w[:, 1] == v[:, 1] .* v[:, 1]
+@test w[:, 2] == v[:, 2] .* v[:, 2]
+@test w[:, 3] == v[:, 3] .* v[:, 3]
 x = copy(v)
 x .= v .* v
 @test x.u == w.u
@@ -153,10 +180,10 @@ w = v .+ 1
 v = DiffEqArray([rand(20), rand(10, 10), rand(3, 3, 3)], 1:3)
 w = v .* v
 @test_broken w isa DiffEqArray # FIXME
-@test w[1] isa Vector
-@test w[1] == v[1] .* v[1]
-@test w[2] == v[2] .* v[2]
-@test w[3] == v[3] .* v[3]
+@test w[:, 1] isa Vector
+@test w[:, 1] == v[:, 1] .* v[:, 1]
+@test w[:, 2] == v[:, 2] .* v[:, 2]
+@test w[:, 3] == v[:, 3] .* v[:, 3]
 x = copy(v)
 x .= v .* v
 @test x.u == w.u
diff --git a/test/downstream/measurements_and_units.jl b/test/downstream/measurements_and_units.jl
index a6a39b60..c3f772f7 100644
--- a/test/downstream/measurements_and_units.jl
+++ b/test/downstream/measurements_and_units.jl
@@ -1,9 +1,9 @@
-using Unitful,MonteCarloMeasurements,OrdinaryDiffEq
+using Unitful, MonteCarloMeasurements, OrdinaryDiffEq
 
 g3 = 9.81u"m/s^2"
-du4 = [10.0±.1,10.0±.1].*u"m/s"
-tspan3 = (0.0,1.0).*u"s"
-f3(du,u,p,t) = [0.0u"m/s^2",-g3]
-u3 = [0.0,0.0].*u"m"
-problem4 = SecondOrderODEProblem(f3,du4,u3,tspan3)
-solve(problem4, Tsit5())
\ No newline at end of file
+du4 = [10.0 ± 0.1, 10.0 ± 0.1] .* u"m/s"
+tspan3 = (0.0, 1.0) .* u"s"
+f3(du, u, p, t) = [0.0u"m/s^2", -g3]
+u3 = [0.0, 0.0] .* u"m"
+problem4 = SecondOrderODEProblem(f3, du4, u3, tspan3)
+@test_broken solve(problem4, Tsit5())
diff --git a/test/downstream/symbol_indexing.jl b/test/downstream/symbol_indexing.jl
index 8110cf7c..29ffdd21 100644
--- a/test/downstream/symbol_indexing.jl
+++ b/test/downstream/symbol_indexing.jl
@@ -1,4 +1,4 @@
-using RecursiveArrayTools, ModelingToolkit, OrdinaryDiffEq, Test
+using RecursiveArrayTools, ModelingToolkit, OrdinaryDiffEq, SymbolicIndexingInterface, Test
 
 include("../testutils.jl")
 
@@ -7,22 +7,23 @@ include("../testutils.jl")
 D = Differential(t)
 @variables RHS(t)
 @named fol_separate = ODESystem([RHS ~ (1 - x) / τ,
-                                    D(x) ~ RHS])
+    D(x) ~ RHS])
 fol_simplified = structural_simplify(fol_separate)
 
 prob = ODEProblem(fol_simplified, [x => 0.0], (0.0, 10.0), [τ => 3.0])
 sol = solve(prob, Tsit5())
 
 sol_new = DiffEqArray(sol.u[1:10],
-                      sol.t[1:10],
-                      sol.prob.f.syms,
-                      sol.prob.f.indepsym,
-                      sol.prob.f.observed,
-                      sol.prob.p)
+    sol.t[1:10],
+    sol.prob.p,
+    sol)
 
 @test sol_new[RHS] ≈ (1 .- sol_new[x]) ./ 3.0
 @test sol_new[t] ≈ sol_new.t
 @test sol_new[t, 1:5] ≈ sol_new.t[1:5]
+@test getp(sol, τ)(sol) == getp(sol_new, τ)(sol_new) == 3.0
+@test_deprecated sol[τ]
+@test_deprecated sol_new[τ]
 
 # Tables interface
 test_tables_interface(sol_new, [:timestamp, Symbol("x(t)")], hcat(sol_new[t], sol_new[x]))
@@ -31,14 +32,14 @@ test_tables_interface(sol_new, [:timestamp, Symbol("x(t)")], hcat(sol_new[t], so
 @variables y(t)
 @parameters α β γ δ
 @named lv = ODESystem([D(x) ~ α * x - β * x * y,
-                          D(y) ~ δ * x * y - γ * x * y])
+    D(y) ~ δ * x * y - γ * x * y])
 
 prob = ODEProblem(lv, [x => 1.0, y => 1.0], (0.0, 10.0),
-                  [α => 1.5, β => 1.0, γ => 3.0, δ => 1.0])
+    [α => 1.5, β => 1.0, γ => 3.0, δ => 1.0])
 sol = solve(prob, Tsit5())
 
 ts = 0:0.5:10
 sol_ts = sol(ts)
 @assert sol_ts isa DiffEqArray
 test_tables_interface(sol_ts, [:timestamp, Symbol("x(t)"), Symbol("y(t)")],
-                      hcat(ts, Array(sol_ts)'))
+    hcat(ts, Array(sol_ts)'))
diff --git a/test/interface_tests.jl b/test/interface_tests.jl
index 3ac5283d..f87bb5dd 100644
--- a/test/interface_tests.jl
+++ b/test/interface_tests.jl
@@ -5,11 +5,11 @@ testva = VectorOfArray([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
 testda = DiffEqArray([[1, 2, 3], [4, 5, 6], [7, 8, 9]], t)
 
 for (i, elem) in enumerate(testva)
-    @test elem == testva[i]
+    @test elem == testva[:, i]
 end
 
 for (i, elem) in enumerate(testda)
-    @test elem == testda[i]
+    @test elem == testda[:, i]
 end
 
 push!(testva, [10, 11, 12])
@@ -43,10 +43,10 @@ push!(testda, [-1, -2, -3, -4])
 @test_throws BoundsError testva[4:5, 5:6]
 @test_throws BoundsError testda[4:5, 5:6]
 
-@test testva[9] == [-1, -2, -3, -4]
-@test testva[end] == [-1, -2, -3, -4]
-@test testda[9] == [-1, -2, -3, -4]
-@test testda[end] == [-1, -2, -3, -4]
+@test testva[:, 9] == [-1, -2, -3, -4]
+@test testva[:, end] == [-1, -2, -3, -4]
+@test testda[:, 9] == [-1, -2, -3, -4]
+@test testda[:, end] == [-1, -2, -3, -4]
 
 # Currently we enforce the general shape, they can just be different lengths, ie we
 # can't do
diff --git a/test/partitions_and_static_arrays.jl b/test/partitions_and_static_arrays.jl
index e26788b4..d5106161 100644
--- a/test/partitions_and_static_arrays.jl
+++ b/test/partitions_and_static_arrays.jl
@@ -1,7 +1,7 @@
 using Test, RecursiveArrayTools, StaticArrays
 
 p = ArrayPartition((zeros(Float32, 2), zeros(SMatrix{2, 2, Int64}, 2),
-                    zeros(SVector{3, Float64}, 2)))
+    zeros(SVector{3, Float64}, 2)))
 @test eltype(p) == Float64
 @test recursive_bottom_eltype(p) == Float64
 @test recursive_unitless_eltype(p) == Float64
diff --git a/test/partitions_test.jl b/test/partitions_test.jl
index cb706638..87bd334a 100644
--- a/test/partitions_test.jl
+++ b/test/partitions_test.jl
@@ -165,7 +165,7 @@ Base.IndexStyle(::MyType) = IndexLinear()
 Base.BroadcastStyle(::Type{<:MyType}) = Broadcast.ArrayStyle{MyType}()
 
 function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{MyType}},
-                      ::Type{T}) where {T}
+    ::Type{T}) where {T}
     similar(find_mt(bc), T)
 end
 
@@ -196,20 +196,20 @@ up = 2 .* ap .+ 1
 @test typeof(ap) == typeof(up)
 
 @testset "ArrayInterface.ismutable(ArrayPartition($a, $b)) == $r" for (a, b, r) in ((1,
-                                                                                     2,
-                                                                                     false),
-                                                                                    ([
-                                                                                         1,
-                                                                                     ],
-                                                                                     2,
-                                                                                     false),
-                                                                                    ([
-                                                                                         1,
-                                                                                     ],
-                                                                                     [
-                                                                                         2,
-                                                                                     ],
-                                                                                     true))
+        2,
+        false),
+    ([
+            1,
+        ],
+        2,
+        false),
+    ([
+            1,
+        ],
+        [
+            2,
+        ],
+        true))
     @test ArrayInterface.ismutable(ArrayPartition(a, b)) == r
 end
 
@@ -231,12 +231,12 @@ begin
     @test_throws MethodError convert(new_type, ArrayPartition(view(b, :), c, c))
 end
 
-
 @testset "Copy and zero with type changing array" begin
     # Motivating use case for this is ArrayPartitions of Arrow arrays which are mmap:ed and change type when copied 
     struct TypeChangingArray{T, N} <: AbstractArray{T, N} end
-    Base.copy(::TypeChangingArray{T, N}) where {T,N} = Array{T,N}(undef, ntuple(_ -> 0, N)) 
-    Base.zero(::TypeChangingArray{T, N}) where {T,N} = zeros(T, ntuple(_ -> 0, N)) 
+    Base.copy(::TypeChangingArray{T, N}) where {T, N} = Array{T, N}(undef,
+        ntuple(_ -> 0, N))
+    Base.zero(::TypeChangingArray{T, N}) where {T, N} = zeros(T, ntuple(_ -> 0, N))
 
     a = ArrayPartition(TypeChangingArray{Int, 2}(), TypeChangingArray{Float32, 2}())
     @test copy(a) == ArrayPartition(zeros(Int, 0, 0), zeros(Float32, 0, 0))
@@ -245,5 +245,5 @@ end
 
 @test !iszero(ArrayPartition([2], [3, 4]))
 @testset "Cartesian indexing" begin
-    @test ArrayPartition([1,2], [3])[1:3,1] == [1, 2, 3]
+    @test ArrayPartition([1, 2], [3])[1:3, 1] == [1, 2, 3]
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 5a45db42..f725cc30 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -26,29 +26,59 @@ end
 
 @time begin
     if GROUP == "Core" || GROUP == "All"
-        @time @safetestset "Utils Tests" begin include("utils_test.jl") end
-        @time @safetestset "Partitions Tests" begin include("partitions_test.jl") end
-        @time @safetestset "VecOfArr Indexing Tests" begin include("basic_indexing.jl") end
-        @time @safetestset "SymbolicIndexingInterface API test" begin include("symbolic_indexing_interface_test.jl") end
-        @time @safetestset "VecOfArr Interface Tests" begin include("interface_tests.jl") end
-        @time @safetestset "Table traits" begin include("tabletraits.jl") end
-        @time @safetestset "StaticArrays Tests" begin include("copy_static_array_test.jl") end
-        @time @safetestset "Linear Algebra Tests" begin include("linalg.jl") end
-        @time @safetestset "Upstream Tests" begin include("upstream.jl") end
-        @time @safetestset "Adjoint Tests" begin include("adjoints.jl") end
-        @time @safetestset "Measurement Tests" begin include("measurements.jl") end
+        @time @safetestset "Utils Tests" begin
+            include("utils_test.jl")
+        end
+        @time @safetestset "Partitions Tests" begin
+            include("partitions_test.jl")
+        end
+        @time @safetestset "VecOfArr Indexing Tests" begin
+            include("basic_indexing.jl")
+        end
+        @time @safetestset "SymbolicIndexingInterface API test" begin
+            include("symbolic_indexing_interface_test.jl")
+        end
+        @time @safetestset "VecOfArr Interface Tests" begin
+            include("interface_tests.jl")
+        end
+        @time @safetestset "Table traits" begin
+            include("tabletraits.jl")
+        end
+        @time @safetestset "StaticArrays Tests" begin
+            include("copy_static_array_test.jl")
+        end
+        @time @safetestset "Linear Algebra Tests" begin
+            include("linalg.jl")
+        end
+        @time @safetestset "Upstream Tests" begin
+            include("upstream.jl")
+        end
+        # @time @safetestset "Adjoint Tests" begin include("adjoints.jl") end
+        @time @safetestset "Measurement Tests" begin
+            include("measurements.jl")
+        end
     end
 
     if !is_APPVEYOR && GROUP == "Downstream"
         activate_downstream_env()
-        @time @safetestset "DiffEqArray Indexing Tests" begin include("downstream/symbol_indexing.jl") end
-        @time @safetestset "Event Tests with ArrayPartition" begin include("downstream/downstream_events.jl") end
-        VERSION >= v"1.9" && @time @safetestset "Measurements and Units" begin include("downstream/measurements_and_units.jl") end
-        @time @safetestset "TrackerExt" begin include("downstream/TrackerExt.jl") end
+        @time @safetestset "DiffEqArray Indexing Tests" begin
+            include("downstream/symbol_indexing.jl")
+        end
+        @time @safetestset "Event Tests with ArrayPartition" begin
+            include("downstream/downstream_events.jl")
+        end
+        VERSION >= v"1.9" && @time @safetestset "Measurements and Units" begin
+            include("downstream/measurements_and_units.jl")
+        end
+        @time @safetestset "TrackerExt" begin
+            include("downstream/TrackerExt.jl")
+        end
     end
 
     if !is_APPVEYOR && GROUP == "GPU"
         activate_gpu_env()
-        @time @safetestset "VectorOfArray GPU" begin include("gpu/vectorofarray_gpu.jl") end
+        @time @safetestset "VectorOfArray GPU" begin
+            include("gpu/vectorofarray_gpu.jl")
+        end
     end
 end
diff --git a/test/symbolic_indexing_interface_test.jl b/test/symbolic_indexing_interface_test.jl
index ad60bbc8..67c3b3b4 100644
--- a/test/symbolic_indexing_interface_test.jl
+++ b/test/symbolic_indexing_interface_test.jl
@@ -1,20 +1,51 @@
-using RecursiveArrayTools, Test, LabelledArrays
+using RecursiveArrayTools, Test, LabelledArrays, SymbolicIndexingInterface
 
 t = 0.0:0.1:1.0
 f(x) = 2x
 f2(x) = 3x
 
-dx = DiffEqArray([[f(x), f2(x)] for x in t], t, [:a, :b], :t)
+dx = DiffEqArray([[f(x), f2(x)] for x in t],
+    t,
+    [1.0, 2.0];
+    variables = [:a, :b],
+    parameters = [:p, :q],
+    independent_variables = [:t])
 @test dx[:t] == t
 @test dx[:a] == [f(x) for x in t]
-@test dx[:b] == [f2(x) for x in t]
-
-dx = DiffEqArray([[f(x), f2(x)] for x in t], t, [:a, :b], [:t])
+@test dx[:a, 2] ≈ f(t[2])
+@test dx[:b, 3] ≈ f2(t[3])
+@test dx[:a, 2:4] ≈ [f(x) for x in t[2:4]]
+@test dx[:b, 4:6] ≈ [f2(x) for x in t[4:6]]
+@test dx[:b] ≈ [f2(x) for x in t]
+@test dx[[:a, :b]] ≈ [[f(x), f2(x)] for x in t]
+@test dx[(:a, :b)] == [(f(x), f2(x)) for x in t]
+@test dx[[:a, :b], 3] ≈ [f(t[3]), f2(t[3])]
+@test dx[[:a, :b], 4:5] ≈ vcat(f.(t[4:5])', f2.(t[4:5])')
+@test getp(dx, [:p, :q])(dx) == [1.0, 2.0]
+@test getp(dx, :p)(dx) == 1.0
+@test getp(dx, :q)(dx) == 2.0
+@test_deprecated dx[:p]
+@test_deprecated dx[[:p, :q]]
 @test dx[:t] == t
-dx = DiffEqArray([[f(x), f2(x)] for x in t], t, [:a, :b])
+
+@test symbolic_container(dx) isa SymbolCache
+@test parameter_values(dx) == [1.0, 2.0]
+@test is_variable.((dx,), [:a, :b, :p, :q, :t]) == [true, true, false, false, false]
+@test variable_index.((dx,), [:a, :b, :p, :q, :t]) == [1, 2, nothing, nothing, nothing]
+@test is_parameter.((dx,), [:a, :b, :p, :q, :t]) == [false, false, true, true, false]
+@test parameter_index.((dx,), [:a, :b, :p, :q, :t]) == [nothing, nothing, 1, 2, nothing]
+@test is_independent_variable.((dx,), [:a, :b, :p, :q, :t]) == [false, false, false, false, true]
+@test variable_symbols(dx) == [:a, :b]
+@test parameter_symbols(dx) == [:p, :q]
+@test independent_variable_symbols(dx) == [:t]
+@test is_time_dependent(dx)
+@test constant_structure(dx)
+
+dx = DiffEqArray([[f(x), f2(x)] for x in t], t; variables = [:a, :b])
 @test_throws Exception dx[nothing] # make sure it isn't storing [nothing] as indepsym
 
 ABC = @SLVector (:a, :b, :c);
 A = ABC(1, 2, 3);
 B = RecursiveArrayTools.DiffEqArray([A, A], [0.0, 2.0]);
 @test getindex(B, :a) == [1, 1]
+
diff --git a/test/testutils.jl b/test/testutils.jl
index 471b65e1..6917b9c7 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -4,7 +4,7 @@ using RecursiveArrayTools: Tables, IteratorInterfaceExtensions
 # Test Tables interface with row access + IteratorInterfaceExtensions for QueryVerse
 # (see https://tables.juliadata.org/stable/#Testing-Tables.jl-Implementations)
 function test_tables_interface(x::AbstractDiffEqArray, names::Vector{Symbol},
-                               values::Matrix)
+    values::Matrix)
     @assert length(names) == size(values, 2)
 
     # AbstractDiffEqArray is a table with row access
diff --git a/test/upstream.jl b/test/upstream.jl
index 467fa8ff..ad79fd94 100644
--- a/test/upstream.jl
+++ b/test/upstream.jl
@@ -33,22 +33,22 @@ nlsolve(mymodel, u0)
 
 function dyn(u, p, t)
     ArrayPartition(ArrayPartition(zeros(1), [0.0]),
-                   ArrayPartition(zeros(1), [0.0]))
+        ArrayPartition(zeros(1), [0.0]))
 end
 
 @test solve(ODEProblem(dyn,
-                       ArrayPartition(ArrayPartition(zeros(1), [-1.0]),
-                                      ArrayPartition(zeros(1), [0.75])),
-                       (0.0, 1.0)), AutoTsit5(Rodas5())).retcode == ReturnCode.Success
+        ArrayPartition(ArrayPartition(zeros(1), [-1.0]),
+            ArrayPartition(zeros(1), [0.75])),
+        (0.0, 1.0)), AutoTsit5(Rodas5())).retcode == ReturnCode.Success
 
 if VERSION < v"1.7"
     @test solve(ODEProblem(dyn,
-                           ArrayPartition(ArrayPartition(zeros(1), [-1.0]),
-                                          ArrayPartition(zeros(1), [0.75])),
-                           (0.0, 1.0)), Rodas5()).retcode == ReturnCode.Success
+            ArrayPartition(ArrayPartition(zeros(1), [-1.0]),
+                ArrayPartition(zeros(1), [0.75])),
+            (0.0, 1.0)), Rodas5()).retcode == ReturnCode.Success
 else
     @test_broken solve(ODEProblem(dyn,
-                                  ArrayPartition(ArrayPartition(zeros(1), [-1.0]),
-                                                 ArrayPartition(zeros(1), [0.75])),
-                                  (0.0, 1.0)), Rodas5()).retcode == ReturnCode.Success
+            ArrayPartition(ArrayPartition(zeros(1), [-1.0]),
+                ArrayPartition(zeros(1), [0.75])),
+            (0.0, 1.0)), Rodas5()).retcode == ReturnCode.Success
 end
diff --git a/test/utils_test.jl b/test/utils_test.jl
index 1b6b6bdd..9769c523 100644
--- a/test/utils_test.jl
+++ b/test/utils_test.jl
@@ -9,7 +9,7 @@ data = convert(Array, randomized)
 ## Test means
 A = [[1 2; 3 4], [1 3; 4 6], [5 6; 7 8]]
 @test recursive_mean(A) ≈ [2.33333333 3.666666666
-                           4.6666666666 6.0]
+    4.6666666666 6.0]
 
 A = zeros(5, 5)
 @test recursive_unitless_eltype(A) == Float64
@@ -67,10 +67,6 @@ function test_recursive_bottom_eltype()
 end
 test_recursive_bottom_eltype()
 
-using RecursiveArrayTools: issymbollike
-@test !issymbollike(1)
-@test issymbollike(:a)
-
 x = zeros(10)
 recursivefill!(x, 1.0)
 @test x == ones(10)
@@ -111,9 +107,9 @@ recursivefill!(x, true)
 
     y_voa = recursivecopy(x_voa)
     recursivefill!(y_voa, true)
-    @test all(y_voa[n] == fill(ones(Vec3), n) for n in 1:4)
+    @test all(y_voa[:, n] == fill(ones(Vec3), n) for n in 1:4)
 
     y_voa = recursivecopy(x_voa)
     recursivefill!(y_voa, ones(Vec3))
-    @test all(y_voa[n] == fill(ones(Vec3), n) for n in 1:4)
+    @test all(y_voa[:, n] == fill(ones(Vec3), n) for n in 1:4)
 end