Different behaviour when plotting one vs many ODE solutions

When taking a slice of the result of solve(ODEProblem), slices with a range (e.g. sol[1:3] return the time series as well, whereas a single element sol[1] only returns the y-values.
I want to be able to plot either one or several solutions in the same line and not sure how to do that.

If I want to plot a number of solutions, I can do:
plot(sol[1:n])

If I want to plot just one, I can do this:
plot(sol.t, sol[1,:])

plot(sol.t, sol[1:n,:]) gives:
ERROR: Expects 8 elements in each col of y, found 3.

and plot(sol[1]) doesn’t plot with the time steps as x values.

Any help appreciated!

Example code below:

using Random, Distributions, OrdinaryDiffEq, Plots
function fun(f, p, t)
    df .= -p[:A] * f
    df
end
rng = MersenneTwister(42)
p = Dict(:A => rand(5,5))
x0 = rand(rng, Uniform(1,10), 5)
tspan=(0,1)
pr = NamedTuple(p)
prob = ODEProblem(fun, x0, tspan, pr)
sol = solve(prob, Tsit5())
# Compare:
# plot(sol[1:3])
#plot(sol.t, sol[1])
# plot(sol.t,sol[1:3,:])
# plot(sol[1])

It plots an interpolation by default.