Recfast splines complain if final redshift is not zero

[jmsull]

If I run the following

𝕡 = CosmoParams()
xstart,xend = -15.0,-1.0
bg = Background(𝕡; x_grid=xstart:0.1:xend)
𝕣 = Bolt.RECFAST(bg=bg, Yp=𝕡.Y_p, OmegaB=𝕡.Ω_b)
ih = IonizationHistory(𝕣, 𝕡, bg)

when I get to the last line I get the following error:

"ERROR: BoundsError: attempt to access 141-element scale(interpolate(OffsetArray(::Array{Float64,1}, 0:142), BSpline(Cubic(Line(OnGrid())))), (-15.0:0.1:-1.0,)) with element type Float64 at index [-0.9980221350275645]"

which is coming from around here since zfinal=0.0 appears hardcoded.

Can be persuaded otherwise but it would be nice if we allowed arbitrary final redshift for the perturbations? I will continue using 0 for now since this is not a priority, but this is something to maybe come back to later.

[xzackli]

@jmsull, I can fix this! I wrote the current RECFAST to replicate the original version as closely as possible, in a way that is not very flexible. This is also related to fixing the differentiability of the RECFAST implementation. The steps would be

1. Crank up the precision, lower the stepsize, and try to find a high-accuracy representation of the solution that the RECFAST code is trying to solve for (it makes various different approximations in different epochs)
2. Rewrite the stepping that is currently done in the dynamic-step-size method that OrdinaryDiffEq.jl uses
3. Make sure that they converge to the same solutions
4. Lower the precision parameters to reach the same level of accuracy as RECFAST, but with the different method

[jmsull]

This all sounds good! One quick thing to mention is that I bumped up Nz to 100000 (this is also the default value in CLASS)