Feat: Support for symbolic indexing of solution objects

[DhairyaLGandhi]

Checklist

• [x] Appropriate tests were added
• [x] Any code changes were done in a way that does not break public API
• [x] All documentation related to code changes were updated
• [ ] The new code follows the contributor guidelines, in particular the SciML Style Guide and COLPRAC.
• [x] Any new documentation only uses public API

Additional context

Allows for the handling of symbolic variables for indexing into solution objects

@parameters σ ρ β
@variables x(t) y(t) z(t) w(t)

eqs = [D(D(x)) ~ σ * (y - x),
    D(y) ~ x * (ρ - z) - y,
    D(z) ~ x * y - β * z,
    w ~ x + y + z]

@mtkbuild sys = ODESystem(eqs, t)

ModelingToolkit.observed(sys)

u0 = [D(x) => 2.0,
    x => 1.0,
    y => 0.0,
    z => 0.0]

p = [σ => 28.0,
    ρ => 10.0,
    β => 8 / 3]

tspan = (0.0, 100.0)
prob = ODEProblem(sys, u0, tspan, p, jac = true)
sol = solve(prob, Tsit5())

sol[sys.x]

[ChrisRackauckas]

Needs a test