Use adaptive k-sampling in full solve command

[hersle]

using SymBoltz, Plots
M = ΛCDM(K = nothing) # flat
pars = SymBoltz.parameters_Planck18(M)
prob = CosmologyProblem(M, pars)

ks = [1.0, 2000.0]
@time sol = solve(prob, ks; adaptiveopts = (rtol = 1e-3,))

τs = sol[M.τ]
Ss = sol(M.ST0, τs, sol.ks)
surface(sol.ks, τs, Ss; ylims = (0.03, 0.08), title = "$(length(sol.ks)) ks")

plot(sol.ks, eachindex(sol.ks); marker = :circle, label = nothing)