Make consistent interface for solution interpolation and indexing

[hersle]

There should be a consistent way for both interpolation and indexing the solution of a CosmologySolution:

• sol(ts, ys) should return the value(s) of the background quantities $y(t)$ on a grid with size length(ts), length(ys). It should interpolate the ODE solution in time.
• sol(ts, ks, ys) should return the value(s) of the perturbation quantities $y(k,t)$ on a grid with size (length(ts), length(ks), length(ys). It should interpolate the ODE solution in time and between wavenumbers.
• sol[it, iy] should return the value(s) of the background quantities $y(t_i)$ at the ODE solver's saved timesteps.
• sol[it, ik, iy] should return the value(s) of the perturbation quantities $y(k_j,t_i)$ at the ODE solver's saved timesteps. It should error if all perturbation ODEs do not have the same time steps.