Age of Universe from Friedmann Equation - How to actually solve the integral? The Friedmann equation for a flat universe can be written as
$$
H(t)=\frac{\dot{a}}{a}=H_0\sqrt{\Omega_{m,0}\cdot a^{-3}+\Omega_{\Lambda,0}}=H(a)
$$
To calculate the age of the universe, many books jump directly to the result. But there should be some sort of integral in between. I assume one can do the following:
$$
t=\int_0^{t_0}\!\mathrm{d}t=\int_0^{1}\!\mathrm{d}a\frac{\mathrm{d}t}{\mathrm{d}a}=\int_0^{1}\!\frac{\mathrm{d}a}{\dot{a}}
$$
with $\dot{a}$ from above expressin for $H$.
But how is this integral solved? Mathematica did something for hours but did not came up with a result. Most books and wikipedia pages skip directly to the result
$$
t_0=\frac{1}{3H_0\sqrt{\Omega_\Lambda}}\log{\frac{1+\sqrt{\Omega_\Lambda}}{1-\sqrt{\Omega_\Lambda}}}
$$
which leads to the well known result of ~13 billion years (depending on the DM density).
Again: But how is the integral solved?
 A: In a flat universe where $\Omega_m+\Omega_{\Lambda}=1$ you can first simplify to
$$H=\frac{\dot{a}}{a}=H_0\sqrt{(1-\Omega_{\Lambda}) a^{-3}+\Omega_{\Lambda}}$$
then you solve for $a$ and get ⁽¹⁾
$$a=\sqrt[3]{\frac{1}{\Omega \Lambda }-1} \ \sinh ^{\frac{2}{3}}\left(\frac{3}{2} H_0 \ t \ \sqrt{\Omega \Lambda }\right)$$
where you solve for $t$ and get
$$t=\frac{2 \sinh ^{-1}\left(\frac{a \sqrt{\frac{a}{\sqrt[3]{\frac{1}{\Omega \Lambda }-1}}}}{\sqrt[3]{\frac{1}{\Omega \Lambda }-1}}\right)}{3 H_0 \sqrt{\Omega \Lambda }}$$
which, when you set $a=1$ for the present scale factor, can be further simplified to
$$t_0=\frac{2 \sinh ^{-1}\left(\frac{1}{\sqrt{\frac{1}{\Omega \Lambda }-1}}\right)}{3 H_0 \sqrt{\Omega \Lambda }}$$
which is equivalent to the expression your source gave.
A: I found that Mathematica could deal with the integral in a minute or so, giving the result:

where I have used the notation $\Omega_{m,0} = m$ and $\Omega_{\Lambda, 0} = n$. Being a conditional expression, it depends upon these parameters being in certain bounds, but let's go ahead and assume that they do. Since the universe is flat we can simplify the result by using the condition $n + m = 1$, and after some steps of algebra we arrive at: $$ H_0 t_0 = \frac{2}{3 \sqrt{n}} \log \frac{1 + \sqrt{n}}{\sqrt{1 - n}} . $$ This is almost the same as the answer you present, except for the denominator of the logarithm, and a factor of 2. Is it possible there is a typo there?
A: for flat universe ($\Omega_m + \Omega_{\Lambda}=1)$
$$H^2 = H_0^2(\Omega_ma^{-3}+\Omega_{\Lambda})$$ 
or $$\frac{\dot{a}^2}{a^2} = H_0(\Omega_ma^{-3}+\Omega_{\Lambda})$$
which becomes
$$\dot{a}^2 = H_0(\Omega_ma^{-1}+\Omega_{\Lambda}a^2)$$
taking square root
$$\dot{a} = H_0\sqrt{\Omega_ma^{-1}+\Omega_{\Lambda}a^2}$$
let us write in the form of
$$\frac{da}{dt} = H_0\sqrt{\frac{\Omega_m}{a}+\Omega_{\Lambda}a^2}$$
$$\frac{da}{\sqrt{\frac{\Omega_m}{a}+\Omega_{\Lambda}a^2}} = H_0dt$$
$$\frac{da}{\sqrt{\Omega_{\Lambda}}\sqrt{(\frac{\Omega_m}{\Omega_{\Lambda}})(\frac{1}{a}) + a^2}}= H_0dt$$
set $$\frac{\Omega_m}{\Omega_{\Lambda}} = L = 0.44927$$
(I took $\Omega_{\Lambda} = 0.69$, $\Omega_m=0.31$)
By taking $a(t_{now})=1$ 
$$\int_0^{a(t_{now})=1}\frac{da}{\sqrt{\frac{L}{a}+a^2}} = \int_0^{t_{universe}} \sqrt{\Omega_{\Lambda}}H_0dt$$
by using wolfram the solution becomes, 
$$\left.\frac{2}{3}log(a^{3/2} + \sqrt{a^2+L})\right|_0^1 = \sqrt{\Omega_{\Lambda}}H_0t_{uni}$$
or you can use https://www.integral-calculator.com to make a numerical calculation of the integral. In any case we have 
$$\int_0^1\frac{da}{\sqrt{\frac{L}{a}+a^2}} = 0.793513$$
Hence,
$$t_{uni} = \frac{0.793513}{0.83066 \times H_0}$$
For $H_0 = 68km/s/Mpc$
$$t_{uni} = 0.9552801 \times H_0^{-1}= 0.9552801 \times 14.39~\text{Gyr} = 13.74 ~\text{Gyr}$$
