La página web que has enlazado es, cuando menos, notable. Contiene toda una cadena de errores que producen un resultado erróneo que los autores interpretan equivocadamente como exacto. Basándome sólo en lo que he visto en esa página, no confiaría en nada más de ese artículo.
Así que, en primer lugar, Wolfram Alpha tiene razón, y te aplaudo por comprobar sus cálculos. Debo señalar que $$ d_{sls} = \frac{c}{H_0\,(1+z_L)}\int_0^{z_L}\left[\Omega_{m,0}(1+z)^3 + \Omega_{\Lambda,0}\right]^{-1/2}\,\text{d}z, $$ así que en tu pregunta has puesto un extra $\sqrt{\Omega_{m,0}}$ en el denominador, que no debería estar ahí.
Veamos dónde van mal las cosas:
El error más flagrante es el paso de (55) a (56): $$ \begin{align} d_{sls} &= \frac{2c}{7H_0(1+z_L)}\left[7\,\Omega_{m,0}^{-1/2} - 2\,\Omega^{\vphantom{1/2}}_{\Lambda,0}\,\Omega_{m,0}^{-3/2} + \mathcal{O}\left((1+z_L)^{-1/2}\right)\right]\\ &\approx \frac{2c}{7H_0(1+z_L)}\left[9\,\Omega_{m,0}^{-1/2} -2\,\Omega_{m,0}^{-3/2}\right] = \frac{2c\,\Omega_{m,0}^{-1/2}}{7H_0(1+z_L)}\left[9 -2\,\Omega_{m,0}^{-1}\right]. \end{align} $$ Así que hay un $\Omega_{m,0}^{-1}$ en lugar de $\Omega_{m,0}^3$ . Como resultado, $d_{sls}/r_s \approx 60$ en lugar de $221$ .
El segundo error es que (55) es incorrecto. La expansión binómica del integrando en $d_{sls}$ sólo es exacta para valores grandes de $z$ pero la integración se realiza en todo el intervalo $[0,z_L]$ y depende fundamentalmente del valor de la función primitiva en $z=0$ . Por eso esta aproximación es tan diferente del cálculo numérico correcto con Wolfram Alpha. Si insertas el valor correcto, deberías obtener $d_{sls}/r_s \approx 149$ (y no $300$ debido a la $\sqrt{\Omega_{m,0}}$ ha añadido usted).
Pero los problemas no acaban ahí. El cálculo del horizonte sonoro $r_s$ también es errónea. La ecuación (47) debería ser $$ r_s = a_L\int_0^{t_L}\frac{c_s}{a(t)}\,\text{d}t. $$ Utilizando (49) y la convención $a_0 = 1$ se puede escribir como $$ r_s = \frac{a_L}{H_0}\int_0^{a_L}\frac{c_s\,\text{d}a}{\sqrt{\Omega_{r,0} + \Omega_{m,0}\,a + \Omega_{k,0}\,a^2 + \Omega_{\Lambda,0}\,a^4}}. $$ A continuación, los autores suponen que la velocidad del sonido $c_s$ permanece más o menos constante y que el término de materia es dominante, de modo que todos los demás términos pueden ignorarse. Esto nos lleva a $$ r_s \approx \frac{c_s\,a_L}{H_0\sqrt{\Omega_{m,0}}}\int_0^{a_L}a^{-1/2}\,\text{d}a = \frac{2\,c_s}{H_0\sqrt{\Omega_{m,0}}}a_L^{3/2} = \frac{2\,c_s}{H_0\sqrt{\Omega_{m,0}}}(1+z_L)^{-3/2}, $$ por lo que existe un factor $3$ errónea en (53). Sin embargo, esto también es una mala aproximación, porque el término de radiación $\Omega_{r,0}$ es bastante significativo en el universo primitivo. Por tanto, una aproximación mejor es $$ \begin{align} r_s &\approx \frac{c_s\,a_L}{H_0}\int_0^{a_L}\frac{\text{d}a}{\sqrt{\Omega_{r,0} + \Omega_{m,0}\,a}}\\ &= \frac{2\,c_s}{H_0\,\Omega_{m,0}}(1+z_L)^{-1}\left[\sqrt{\Omega_{r,0} + \Omega_{m,0}\,(1+z_L)^{-1}} - \sqrt{\Omega_{r,0}}\,\right]. \end{align} $$ Utilizando el valor $\Omega_{r,0}\approx 9\times 10^{-5}$ Esto conduce a $d_{sls}/r_s \approx 86$ . Un cálculo completo con $c_s$ que no derivaré aquí, da $d_{sls}/r_s \approx 96$ .
Por último, una aproximación de primer orden a la posición del primer pico acústico no viene dada por $d_{sls}/r_s$ sino por $\pi d_{sls}/r_s$ y encontramos $\pi 96 \approx 300$ . Como ves, esto todavía no está muy cerca del valor real de $\approx 220$ . La razón es que los picos acústicos de bajo orden están muy influidos por la propagación de las fluctuaciones de densidad y por los efectos gravitatorios de corrimiento al rojo en el universo primitivo. El tratamiento completo de estos es muy complejo, y sólo para especialistas (que yo no soy, por cierto).
Creo que nunca he visto un sitio web que parezca profesional, pero que esté tan mal en tantos niveles. Mi consejo sería encontrar un curso mejor. Por ejemplo, el libro de texto de cosmología de Weinberg tiene todo lo que necesitas, pero no es una lectura fácil.