We donâ€™t have an initial condition for the ODE for the optical depth \(\tau\) if we solve it from the early Universe till today. However we can just pick any initial value, solve the ODE and then normalize it by subtracting off the value at \(x=0\) to get a function that has \( \tau(0) = 0 \). This works fine, but there is one issue to be aware of. The optical depth we subtract off doing this normalization can be huge, lets say, \(10^{10} \). The resulting function \(\tau = \tau_{\rm ODE} - 10^{10}\) will then have \(\tau(0) \equiv 0\) as we want, but it will only have an (absolute) numerical accuracy of just \(\sim 10^{-6}\) in this case (double precision means a relative precision of \(10^{-16}\)). Thus close to \(x=0\) where \(\tau\) is small we might get some numerical noise popping up. This might not be seen in the function itself, but it will be seen in the derivatives which is more sensitive to this.

**Solutions:**

Fill an array with \(d\tau(x)/dx\) just as you do in the ODE solver using the formula for this. Make a spline of this and use this spline to get \(d\tau(x)/dx\) and its derivative to get \(d^2\tau(x)/dx^2\). This usually works very well. Another option is to just solve the ODE from today and back in time (I just updated the ODEsolver on GitHub to allow this) and then its no issue.