diff --git a/src/observables/fourier.jl b/src/observables/fourier.jl
index 3d2d3191..1a065df8 100644
--- a/src/observables/fourier.jl
+++ b/src/observables/fourier.jl
@@ -80,8 +80,14 @@ end
 Solve the problem `prob` with exact wavenumber(s) `k`, and then compute the power spectrum with the solution `sol`.
 """
 function spectrum_matter(prob::CosmologyProblem, k, τ = nothing; species = [:c, :b, :h], kwargs...)
-    sol = solve(prob, k; kwargs...) # TODO: just save endpoints
-    τ = isnothing(τ) ? sol[prob.M.τ][end] : τ
+    # save only the necessary time(s)
+    if isnothing(τ)
+        ptextraopts = (save_everystep = false, save_start = false, save_end = true)
+    else
+        ptextraopts = (saveat = [τ],)
+    end
+    sol = solve(prob, k; ptextraopts, kwargs...)
+    τ = isnothing(τ) ? sol.bg.t[end] : τ
     return spectrum_matter(sol, k, τ; species)
 end
 
diff --git a/src/solve.jl b/src/solve.jl
index 0fba9d04..b2942034 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -267,8 +267,8 @@ end
 """
     function solve(
         prob::CosmologyProblem, ks::Union{Nothing, AbstractArray} = nothing;
-        bgopts = (alg = DEFAULT_BGALG, reltol = 1e-9, abstol = 1e-9),
-        ptopts = (alg = DEFAULT_PTALG, reltol = 1e-8, abstol = 1e-8),
+        bgopts = (alg = DEFAULT_BGALG, reltol = 1e-9, abstol = 1e-9), bgextraopts = (),
+        ptopts = (alg = DEFAULT_PTALG, reltol = 1e-8, abstol = 1e-8), ptextraopts = (),
         shootopts = (alg = DEFAULT_SHOOTALG, abstol = 1e-5),
         thread = true, verbose = false, kwargs...
     )
@@ -280,15 +280,15 @@ If `threads`, integration over independent perturbation modes are parallellized.
 """
 function solve(
     prob::CosmologyProblem, ks::Union{Nothing, AbstractArray} = nothing;
-    bgopts = (alg = DEFAULT_BGALG, reltol = 1e-9, abstol = 1e-9),
-    ptopts = (alg = DEFAULT_PTALG, reltol = 1e-8, abstol = 1e-8),
+    bgopts = (alg = DEFAULT_BGALG, reltol = 1e-9, abstol = 1e-9), bgextraopts = (),
+    ptopts = (alg = DEFAULT_PTALG, reltol = 1e-8, abstol = 1e-8), ptextraopts = (),
     shootopts = (alg = DEFAULT_SHOOTALG, abstol = 1e-5),
     thread = true, verbose = false, kwargs...
 )
     if !isempty(prob.shoot)
-        bgsol = solvebg(prob.bg, collect(prob.shoot), prob.conditions; shootopts, verbose, bgopts..., kwargs...)
+        bgsol = solvebg(prob.bg, collect(prob.shoot), prob.conditions; shootopts, verbose, bgopts..., bgextraopts..., kwargs...)
     else
-        bgsol = solvebg(prob.bg; verbose, bgopts..., kwargs...)
+        bgsol = solvebg(prob.bg; verbose, bgopts..., bgextraopts..., kwargs...)
     end
 
     if isnothing(ks) || isempty(ks)
@@ -296,7 +296,7 @@ function solve(
         ptsol = nothing
     else
         ks = k_dimensionless.(ks, Ref(bgsol))
-        ptsol = solvept(prob.pt, bgsol, ks, prob.bgspline; thread, verbose, ptopts..., kwargs...)
+        ptsol = solvept(prob.pt, bgsol, ks, prob.bgspline; thread, verbose, ptopts..., ptextraopts..., kwargs...)
     end
 
     return CosmologySolution(prob, bgsol, ks, ptsol)