diff --git a/docs/pages.jl b/docs/pages.jl
index 5d9ec136..16338f15 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -2,6 +2,7 @@
 
 pages = [
     "Home" => "index.md",
+    "AbstractVectorOfArrayInterface.md",
     "array_types.md",
     "recursive_array_functions.md"
 ]
diff --git a/docs/src/AbstractVectorOfArrayInterface.md b/docs/src/AbstractVectorOfArrayInterface.md
new file mode 100644
index 00000000..907cc511
--- /dev/null
+++ b/docs/src/AbstractVectorOfArrayInterface.md
@@ -0,0 +1,6 @@
+# The AbstractVectorOfArray and AbstractDiffEqArray Interfaces
+
+```@doc
+AbstractVectorOfArray
+AbstractDiffEqArray
+```
diff --git a/src/RecursiveArrayTools.jl b/src/RecursiveArrayTools.jl
index 9e251ac0..a616c592 100644
--- a/src/RecursiveArrayTools.jl
+++ b/src/RecursiveArrayTools.jl
@@ -13,7 +13,120 @@ import Adapt
 
 import Tables, IteratorInterfaceExtensions
 
+"""
+    AbstractVectorOfArray{T, N, A}
+
+An AbstractVectorOfArray is an object which represents arrays of arrays,
+and arbitrary recursive nesting of arrays, as a single array-like object.
+Thus a canonical example of an AbstractVectorOfArray is something of the
+form `VectorOfArray([[1,2],[3,4]])`, which "acts" like the matrix [1 3; 2 4]
+where the data is stored and accessed in a column-ordered fashion (as is typical
+in Julia), but the actual matrix is never constructed and instead lazily represented
+through the type.
+
+An AbstractVectorOfArray subtype should match the following behaviors.
+
+!!! note
+
+    In 2023 the linear indexing `A[i]`` was deprecated. It previously had the behavior that
+    `A[i] = A.u[i]`. However, this is incompatible with standard `AbstractArray` interfaces,
+    Since if `A = VectorOfArray([[1,2],[3,4]])` and `A` is supposed to act like `[1 3; 2 4]`,
+    then there is a difference `A[1] = [1,2]` for the VectorOfArray while `A[1] = 1` for the
+    matrix. This causes many issues if `AbstractVectorOfArray <: AbstractArray`. Thus we
+    plan in 2026 to complete the deprecation and thus have a breaking update where `A[i]`
+    matches the linear indexing of an `AbstractArray`, and then making
+    `AbstractVectorOfArray <: AbstractArray`. Until then, `AbstractVectorOfArray` due to
+    this interface break but manaully implements an AbstractArray-like interface for
+    future compatability.
+
+## Fields
+
+An AbstractVectorOfArray has the following fields:
+
+* `u` which holds the Vector of values at each timestep
+
+## Array Interface
+
+The general operations are as follows. Use
+
+```julia
+A.u[j]
+```
+
+to access the `j`th array. For multidimensional systems, this
+will address first by component and lastly by time, and thus
+
+```julia
+A[i, j]
+```
+
+will be the `i`th component at array `j`. Hence, `A[j][i] == A[i, j]`. This is done
+because Julia is column-major, so the leading dimension should be contiguous in memory.
+If the independent variables had shape (for example, was a matrix), then `i` is the
+linear index. We can also access solutions with shape:
+
+```julia
+A[i, k, j]
+```
+
+gives the `[i,k]` component of the system at array `j`. The colon operator is
+supported, meaning that
+
+```julia
+A[i, :]
+```
+
+gives the timeseries for the `i`th component.
+
+## Using the AbstractArray Interface
+
+The `AbstractArray` interface can be directly used. For example, for a vector
+system of variables `A[i,j]` is a matrix with rows being the variables and
+columns being the timepoints. Operations like `A'` will
+transpose the solution type. Functionality written for `AbstractArray`s can
+directly use this. For example, the Base `cov` function computes correlations
+amongst columns, and thus:
+
+```julia
+cov(A)
+```
+
+computes the correlation of the system state in time, whereas
+
+```julia
+cov(A, 2)
+```
+
+computes the correlation between the variables. Similarly, `mean(A,2)` is the
+mean of the variable in time, and `var(A,2)` is the variance. Other statistical
+functions and packages which work on `AbstractArray` types will work on the
+solution type.
+
+## Conversions
+
+At anytime, a true `Array` can be created using `Array(A)`, or more generally `stack(A)`
+to make the array type match the internal array type (for example, if `A` is an array
+of GPU arrays, `stack(A)` will be a GPU array).
+"""
 abstract type AbstractVectorOfArray{T, N, A} end
+
+"""
+    AbstractDiffEqArray{T, N, A} <: AbstractVectorOfArray{T, N, A}
+
+An AbstractVectorOfArray object which has extra information of a time array `A.t`
+in order to specify a time series. A canonical AbstractDiffEqArray is for example
+the pairing `DiffEqArray([[1,2],[3,4]],[1.0,2.0])` which means that at time 1.0
+the values were `[1,2]` and at time 2.0 the values were `[3,4]`.
+
+An AbstractDiffEqArray has all of the same behaviors as an AbstractVectorOfArray with the
+additional properties:
+
+## Fields
+
+An AbstractDiffEqArray adds the following fields:
+
+* `t` which holds the times of each timestep.
+"""
 abstract type AbstractDiffEqArray{T, N, A} <: AbstractVectorOfArray{T, N, A} end
 
 include("utils.jl")
diff --git a/src/vector_of_array.jl b/src/vector_of_array.jl
index 7e38fb89..c05a35c2 100644
--- a/src/vector_of_array.jl
+++ b/src/vector_of_array.jl
@@ -388,29 +388,32 @@ end
 @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(
+function Base.firstindex(VA::AbstractVectorOfArray{T,N,A}) where {T,N,A}
+    N > 1 && Base.depwarn(
         "Linear indexing of `AbstractVectorOfArray` is deprecated. Change `A[i]` to `A.u[i]` ",
         :firstindex)
     return firstindex(VA.u)
 end
 
-function Base.lastindex(VA::AbstractVectorOfArray)
-    Base.depwarn(
+function Base.lastindex(VA::AbstractVectorOfArray{T,N,A}) where {T,N,A}
+     N > 1 && 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
+Base.getindex(A::AbstractVectorOfArray, I::Int) = A.u[I]
+Base.getindex(A::AbstractVectorOfArray, I::AbstractArray{Int}) = A.u[I]
+Base.getindex(A::AbstractDiffEqArray, I::Int) = A.u[I]
+Base.getindex(A::AbstractDiffEqArray, I::AbstractArray{Int}) = A.u[I]
 
-@deprecate Base.getindex(A::AbstractVectorOfArray, I::AbstractArray{Int}) Base.getindex(
-    A, :, I) false
+@deprecate Base.getindex(VA::AbstractVectorOfArray{T,N,A}, I::Int) where {T,N,A<:Union{AbstractArray, AbstractVectorOfArray}} VA.u[I] false
 
-@deprecate Base.getindex(A::AbstractDiffEqArray, I::AbstractArray{Int}) Base.getindex(
-    A, :, I) false
+@deprecate Base.getindex(VA::AbstractVectorOfArray{T,N,A}, I::AbstractArray{Int}) where {T,N,A<:Union{AbstractArray, AbstractVectorOfArray}} VA.u[I] false
 
-@deprecate Base.getindex(A::AbstractDiffEqArray, i::Int) Base.getindex(A, :, i) false
+@deprecate Base.getindex(VA::AbstractDiffEqArray{T,N,A}, I::AbstractArray{Int}) where {T,N,A<:Union{AbstractArray, AbstractVectorOfArray}} VA.u[I] false
+
+@deprecate Base.getindex(VA::AbstractDiffEqArray{T,N,A}, i::Int) where {T,N,A<:Union{AbstractArray, AbstractVectorOfArray}} VA.u[i] false
 
 __parameterless_type(T) = Base.typename(T).wrapper
 
@@ -520,22 +523,25 @@ Base.@propagate_inbounds function Base.setindex!(VA::AbstractVectorOfArray{T, N}
     VA.u[I] = v
 end
 
-@deprecate Base.setindex!(VA::AbstractVectorOfArray, v, I::Int) Base.setindex!(VA, v, :, I) false
+Base.@propagate_inbounds Base.setindex!(VA::AbstractVectorOfArray, v, I::Int) = Base.setindex!(VA.u, v, I)
+@deprecate Base.setindex!(VA::AbstractVectorOfArray{T,N,A}, v, I::Int) where {T,N,A<:Union{AbstractArray, AbstractVectorOfArray}} Base.setindex!(VA.u, v, I) false
 
 Base.@propagate_inbounds function Base.setindex!(VA::AbstractVectorOfArray{T, N}, v,
         ::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 Base.setindex!(VA::AbstractVectorOfArray, v, I::Colon) = Base.setindex!(VA.u, v, I)
+@deprecate Base.setindex!(VA::AbstractVectorOfArray{T,N,A}, v, I::Colon)  where {T,N,A<:Union{AbstractArray, AbstractVectorOfArray}} Base.setindex!(
+    VA.u, v, I) false
 
 Base.@propagate_inbounds function Base.setindex!(VA::AbstractVectorOfArray{T, N}, v,
         ::Colon, I::AbstractArray{Int}) where {T, N}
     VA.u[I] = v
 end
 
-@deprecate Base.setindex!(VA::AbstractVectorOfArray, v, I::AbstractArray{Int}) Base.setindex!(
+Base.@propagate_inbounds Base.setindex!(VA::AbstractVectorOfArray, v, I::AbstractArray{Int}) = Base.setindex!(VA.u, v, I)
+@deprecate Base.setindex!(VA::AbstractVectorOfArray{T,N,A}, v, I::AbstractArray{Int}) where {T,N,A<:Union{AbstractArray, AbstractVectorOfArray}} Base.setindex!(
     VA, v, :, I) false
 
 Base.@propagate_inbounds function Base.setindex!(