diff --git a/src/solve.jl b/src/solve.jl
index 6f3d0c8d..62fc5001 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -145,9 +145,11 @@ function solve(M::CosmologyModel, pars; aini = 1e-7, solver = Rodas5P(), reltol
         check_solution(th_sol.retcode)
 
         # Offset optical depth, so it's 0 today
-        idx_τ = variable_index(M.th, M.b.rec.τ)
-        for i in 1:length(th_sol.u)
-            th_sol.u[i][idx_τ] -= th_sol.u[end][idx_τ]
+        if have(M.sys, :b)
+            idx_τ = variable_index(M.th, M.b.rec.τ)
+            for i in 1:length(th_sol.u)
+                th_sol.u[i][idx_τ] -= th_sol.u[end][idx_τ]
+            end
         end
     else
         th_sol = bg_sol
diff --git a/src/utils.jl b/src/utils.jl
index 80b2db87..ebe07b7c 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -71,6 +71,7 @@ function taylor(sys::ODESystem, ϵ, orders)
     return ODESystem(eqs, t, vars, pars; initialization_eqs=ieqs, defaults=defs, guesses=guesses, name=sys.name, description=sys.description)
 end
 
+have(sys, s) = s in nameof.(ModelingToolkit.get_systems(sys))
 have(s) = !isnothing(s) # shorthand for checking if we have a given species
 
 function spline(y, x)