The equation using Hubble velocity along the line of sight to do the integral:

τu0=πe2λ0mecf12nHIπbexp[(uu0b)2]1Hdv, where dv=Hdr is the differential in Hubble velocity along the LOS

For the second form of the equation, using the gas velocity Hr+upec to do the integral, then

dr=|dudr|1du=|1H+upec|du

Where

du=|Hdr+Δupec|

is the differential of the gas velocity along the LOS. Note that an absolute value around du is needed to avoid negative values.

τu0=πe2λ0mecf12nHIπbexp[(uu0b)2]|1H+upec|du

Numerically the gradient and the difference are evaluated:

upec=ui+1pecui1pec2dr Δupec=ui+1pecui1pec2

The comparison of the Flux from the two calculations is below.

The fractional differences in the Flux are 1011!!!

Actually, the term

Hdr+ΔupecH+upec=Hdr+ΔupecH+Δupecdr=drH+ΔupecrH+Δupecdr=dr

Which of course should hold from the Jacobian of the transformation. Then this becomes the original form of the equation just changing dr=H1dv

The issue that I had before was that I was missing the term Δupec in du.