BREAKING: Add ragged arrays, interpolation interface, and plotting delegation

## Ragged Arrays (RecursiveArrayToolsRaggedArrays sublibrary)

- Add AbstractRaggedVectorOfArray and AbstractRaggedDiffEqArray abstract types
- Add RaggedVectorOfArray and RaggedDiffEqArray concrete types that preserve
  true ragged structure without zero-padding
- Column-major linear indexing, no-padding A[:, i], symbolic indexing,
  broadcasting, and callable interpolation interface
- Isolated in lib/ sublibrary to avoid invalidations

## Interpolation Interface on DiffEqArray

- Add `interp` and `dense` fields to DiffEqArray (new type parameter I)
- Callable syntax: da(t), da(t; idxs=1), da(t, Val{1})
- All constructors accept `interp=nothing, dense=false` kwargs

## Plotting Delegation

- Full-featured AbstractDiffEqArray plot recipe with:
  - idxs variable selection (integer, array, tuple for phase plots, symbolic)
  - Dense interpolation via callable interface
  - denseplot, plotdensity, tspan, plotat keyword arguments
- Export plotting helpers for SciMLBase delegation:
  DEFAULT_PLOT_FUNC, plottable_indices, plot_indices, getindepsym_defaultt,
  interpret_vars, add_labels!, diffeq_to_arrays, solplot_vecs_and_labels

## Bug Fix

- Fix ragged VectorOfArray .= zero broadcast (DimensionMismatch)
  by using dest.u[i] instead of dest[:, i] in copyto!

## Version

- Bump to v4.0.0

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: ChrisRackauckas

Project.toml

@@ -1,6 +1,6 @@
 name = "RecursiveArrayTools"
 uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
-version = "3.48.0"
+version = "4.0.0"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
 
 [deps]
@@ -11,6 +11,8 @@ GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SnoopCompileCore = "e2b509da-e806-4183-be48-004708413034"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 
@@ -62,6 +64,7 @@ RecipesBase = "1.3.4"
 ReverseDiff = "1.15"
 SafeTestsets = "0.1"
 SciMLBase = "2.103"
+SnoopCompileCore = "3.1.1"
 SparseArrays = "1.10"
 StaticArrays = "1.6"
 StaticArraysCore = "1.4.2"

docs/pages.jl

@@ -5,5 +5,6 @@ pages = [
     "breaking_changes_v4.md",
     "AbstractVectorOfArrayInterface.md",
     "array_types.md",
+    "plotting.md",
     "recursive_array_functions.md",
 ]

docs/src/breaking_changes_v4.md

@@ -82,6 +82,40 @@ not indices into `A.u`.
 | `map(f, A)` (map over columns) | `map(f, A.u)` |
 | `A == vec_of_vecs` | `A.u == vec_of_vecs` |
 
+## Interpolation Interface on DiffEqArray
+
+`DiffEqArray` now has `interp` and `dense` fields for interpolation support:
+
+```julia
+# Create with interpolation
+da = DiffEqArray(u, t, p, sys; interp = my_interp, dense = true)
+
+# Callable syntax
+da(0.5)                    # interpolate at t=0.5
+da(0.5, Val{1})            # first derivative at t=0.5
+da([0.1, 0.5, 0.9])       # interpolate at multiple times
+da(0.5; idxs = 1)          # interpolate single component
+da(0.5; idxs = [1, 2])    # interpolate subset of components
+```
+
+The interpolation object must be callable as `interp(t, idxs, deriv, p, continuity)`,
+matching the protocol used by SciMLBase's `LinearInterpolation`, `HermiteInterpolation`,
+and `ConstantInterpolation`.
+
+When `dense = true` and `interp` is provided, `plot(da)` automatically generates
+dense interpolated output instead of plotting only the saved time points.
+
+## Ragged Arrays Sublibrary
+
+`RaggedVectorOfArray` and `RaggedDiffEqArray` are available via:
+
+```julia
+using RecursiveArrayToolsRaggedArrays
+```
+
+These types preserve the true ragged structure without zero-padding, and do **not**
+subtype `AbstractArray`. See the `RecursiveArrayToolsRaggedArrays` subpackage for details.
+
 ## Zygote Compatibility
 
 Some Zygote adjoint rules need updating for the new `AbstractArray` subtyping.

docs/src/plotting.md

@@ -0,0 +1,154 @@
+# Plotting
+
+`AbstractVectorOfArray` and `AbstractDiffEqArray` types include
+[Plots.jl](https://github.com/JuliaPlots/Plots.jl) recipes so they can be
+visualised directly with `plot(A)`.
+
+## VectorOfArray
+
+A `VectorOfArray` plots as a matrix where each inner array is a column:
+
+```julia
+using RecursiveArrayTools, Plots
+A = VectorOfArray([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+plot(A)   # 3 series (components) vs column index
+```
+
+## DiffEqArray
+
+A `DiffEqArray` plots as component time series against `A.t`:
+
+```julia
+u = [[sin(t), cos(t)] for t in 0:0.1:2pi]
+t = collect(0:0.1:2pi)
+A = DiffEqArray(u, t)
+plot(A)   # plots sin and cos vs t
+```
+
+If the `DiffEqArray` carries a symbolic system (via `variables` and
+`independent_variables` keyword arguments), the axis labels and series names
+are set automatically from the symbol names.
+
+## Dense (Interpolated) Plotting
+
+When a `DiffEqArray` has an interpolation object in its `interp` field and
+`dense = true`, calling `plot(A)` generates a smooth curve by evaluating the
+interpolation at many points rather than connecting only the saved time steps.
+
+```julia
+using SciMLBase: LinearInterpolation
+
+u = [[1.0, 0.0], [0.0, 1.0], [-1.0, 0.0], [0.0, -1.0]]
+t = [0.0, 1.0, 2.0, 3.0]
+interp = LinearInterpolation(t, u)
+A = DiffEqArray(u, t; interp = interp, dense = true)
+plot(A)   # smooth interpolated curve with 1000+ points
+```
+
+### Plot Recipe Keyword Arguments
+
+The `AbstractDiffEqArray` recipe accepts the following keyword arguments,
+which can be passed directly to `plot`:
+
+| Keyword | Type | Default | Description |
+|---------|------|---------|-------------|
+| `denseplot` | `Bool` | `A.dense && A.interp !== nothing` | Use dense interpolation for smooth curves. Set `false` to show only saved points. |
+| `plotdensity` | `Int` | `max(1000, 10 * length(A.u))` | Number of evenly-spaced points to evaluate when `denseplot = true`. |
+| `tspan` | `Tuple` or `nothing` | `nothing` | Restrict the time window. E.g. `tspan = (0.0, 5.0)`. |
+| `idxs` | varies | `nothing` | Select which components to plot (see below). |
+
+Example:
+
+```julia
+plot(A; denseplot = true, plotdensity = 5000, tspan = (0.0, 2.0))
+```
+
+## Selecting Variables with `idxs`
+
+The `idxs` keyword controls which variables appear in the plot. It supports
+several formats:
+
+### Single component
+
+```julia
+plot(A; idxs = 1)      # plot component 1 vs time
+plot(A; idxs = 2)      # plot component 2 vs time
+```
+
+### Multiple components
+
+```julia
+plot(A; idxs = [1, 3, 5])   # plot components 1, 3, 5 vs time
+```
+
+### Phase-space plots (component vs component)
+
+Use a tuple where index `0` represents the independent variable (time):
+
+```julia
+plot(A; idxs = (1, 2))      # component 1 vs component 2
+plot(A; idxs = (0, 1))      # time vs component 1 (same as default)
+plot(A; idxs = (1, 2, 3))   # 3D plot of components 1, 2, 3
+```
+
+### Symbolic indexing
+
+When the `DiffEqArray` carries a symbolic system, variables can be referenced
+by symbol:
+
+```julia
+A = DiffEqArray(u, t; variables = [:x, :y], independent_variables = [:t])
+plot(A; idxs = :x)           # plot x vs time
+plot(A; idxs = [:x, :y])     # plot both
+plot(A; idxs = (:x, :y))     # phase plot of x vs y
+```
+
+### Custom transformations
+
+A function can be applied to the plotted values:
+
+```julia
+plot(A; idxs = (norm, 0, 1, 2))   # plot norm(u1, u2) vs time
+```
+
+The tuple format is `(f, xvar, yvar)` or `(f, xvar, yvar, zvar)` where `f`
+is applied element-wise.
+
+## Callable Interface
+
+Any `AbstractDiffEqArray` with an `interp` field supports callable syntax for
+interpolation, independent of plotting:
+
+```julia
+A(0.5)                    # interpolate all components at t=0.5
+A(0.5; idxs = 1)          # interpolate component 1 at t=0.5
+A([0.1, 0.5, 0.9])       # interpolate at multiple times (returns DiffEqArray)
+A(0.5, Val{1})            # first derivative at t=0.5
+A(0.5; continuity = :right)  # right-continuity at discontinuities
+```
+
+The interpolation object must be callable as
+`interp(t, idxs, deriv, p, continuity)`, matching the protocol used by
+SciMLBase's `LinearInterpolation`, `HermiteInterpolation`, and
+`ConstantInterpolation`.
+
+When no interpolation is available (`interp === nothing`), calling `A(t)`
+throws an error.
+
+## ODE Solution Plotting
+
+ODE solutions from DifferentialEquations.jl are subtypes of
+`AbstractDiffEqArray` and inherit all of the above functionality, plus
+additional features:
+
+- **Automatic dense plotting**: `denseplot` defaults to `true` when the solver
+  provides dense output.
+- **Analytic solution overlay**: `plot(sol; plot_analytic = true)` overlays the
+  exact solution if the problem defines one.
+- **Discrete variables**: Time-varying parameters are plotted as step functions
+  with dashed lines and markers.
+- **Symbolic indexing**: Full symbolic indexing through ModelingToolkit is
+  supported, including observed (derived) variables.
+
+These advanced features are defined in SciMLBase and activate automatically
+when plotting solution objects.

lib/RecursiveArrayToolsRaggedArrays/Project.toml

@@ -0,0 +1,19 @@
+name = "RecursiveArrayToolsRaggedArrays"
+uuid = "c384ba91-639a-44ca-823a-e1d3691ab84a"
+version = "1.0.0"
+
+[deps]
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
+
+[compat]
+RecursiveArrayTools = "3.48"
+SymbolicIndexingInterface = "0.3.35"
+julia = "1.10"
+
+[extras]
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["SymbolicIndexingInterface", "Test"]