Let's consider that $\vec{r}_0$ is a position of the Earth and $\vec{x}$ is a position of a certain point in the universe. So $\vec{r} = \vec{x} - \vec{r}_0$ is a position relative to the Earth.
According to the Hubble law a time derivative $\dot{\vec{R}}$ of any vector $\vec{R}$ is $H\vec{R}$. Thus
\[ \dot{\vec{r}}= - \dot{\vec{r}}_0 + \dot{\vec{x}} = H \vec{r} \]
and the Hubble constant $H$ is independent of the observer's position.
Remark: if the mass of the system increases from $m$ to $m + \Delta m$ without any changing the relative mass distribution, the gravitational potential energy increases by
\[ \Delta E_G = \int\limits_0^{\Delta m} \varphi \, dm,\]
where $\varphi$ is a gravitational potential at the point of mass $dm$. Also $\varphi$ is chosen so that $\varphi=0$ on infinity.
Gauss's law for the gravitational field has the form \[\int_\Gamma \vec{g} \cdot d\vec{S} = -4\pi Gm,\] where $m$ is the mass enclosed within the closed surface $\Gamma$. Let's choose a sphere with radius $r$ and center coinciding with the center of the universe. So
\[ - 4\pi r^2 g = -4\pi G m\frac{r^3}{R^3} \quad \Rightarrow \quad g = GM \frac{r}{R^3}.\]
Gravitational potential at the surface of the universe $\varphi(R) = -GM/R$, so
\[ \varphi(r) = -\frac{GM}{R} - \frac{GM}{R^3}\int\limits_r^R r \, dr = - \frac{3GM}{2R} + \frac{GMr^2}{2R^3}.\]
Now we can calculate the change $\Delta E_G$ in the gravitational potential energy as $M \to M + \Delta M$.
\[\Delta E_G = \int\limits_0^{R} \left( - \frac{3GM}{2R} + \frac{GMr^2}{2R^3} \right) \frac{4\pi r^2 dr}{\frac{4}{3}\pi R^3} \Delta M = -\frac{3GM\Delta M}{2R}+\frac{3GM\Delta M}{10R} = -\frac{6 GM \Delta M}{5R}.\]
Thus
\[E_G = -\frac{6G}{5R} \int\limits_0^{M} m \, dm = -\frac{3GM^2}{5R}.\]
We can calculate $E_K$ by integrating:
\[E_K = \int\limits_0^{R} \frac{(Hr)^2}{2} \cdot M \frac{4\pi r^2 \, dr}{\frac{4}{3} \pi R^3} = \frac{3H^2 M}{2R^3} \int\limits_0^R r^4 \, dr = \frac{3H^2MR^2}{10} \]
Total energy of the universe
\[ E = -\frac{3GM^2}{5R} + \frac{3H^2MR^2}{10},\]
where $M=4\pi \rho R^3/3,$ so
\[ E = -\frac{16 \pi^2}{15}\rho^2 GR^5 + \frac{2\pi}{5} H^2 \rho R^5. \]
If $E<0$ the universe has maximum size when $E = E_G$, so
\[ E = \frac{2\pi}{5} \rho R^5 \left[ - \frac{8 \pi}{3} \rho G + H^2 \right] < 0 \quad \Rightarrow \quad H^2 < \frac{8 \pi}{3} \rho G \]
Maximum radius is reached when $E=E_G$, so
\[-\frac{2}{15}Mc^2 = -\frac{3GM^2}{5R_\text{max}} \quad \Rightarrow \quad R_\text{max} = \frac{9GM}{2c^2}\]
Remark: $$ \int\limits_{0}^{1} \frac{\sqrt{x}}{\sqrt{1-x}} d x=(y=\sqrt{1-x})=2 \int\limits_{0}^{1} \sqrt{1-y^{2}} d y $$
As we discussed earlier, the total energy of the universe
\[ E = -\frac{3GM^2}{5R} + \frac{3H^2MR^2}{10},\]
where $\dot{R} = HR$. Thus,
\[ \dot{R} = \pm \sqrt{\frac{2GM}{R} + \frac{10E}{3M}} = \pm\sqrt{2GM \left(\frac{1}{R} - \frac{1}{R_\text{max}} \right) }.\]
The lifetime consists of the expansion time and the compression time, which are equal. Let's find the expansion time:
\[T_e = \int\limits_0^{R_\text{max}} \frac{dR}{\sqrt{2GM\left(\frac{1}{R} - \frac{1}{R_\text{max}} \right)}} = \sqrt{\frac{R_\text{max}^3}{2GM}} \int\limits_0^1 \frac{\sqrt{\frac{R}{R_\text{max}}} \, d\left(\frac{R}{R_\text{max}} \right)}{\sqrt{1 - \frac{R}{R_\text{max}}}}.\]
To calculate the integral mentioned in the remark, let's use the substitution $y = \sin \phi$, $dy = \cos \phi \, d \phi$:
\[2 \int\limits_0^1 \sqrt{1-y^2}dy = 2 \int\limits_0^{\pi/2}\cos^2 \phi \, d\phi = \int\limits_0^{\pi/2} \left( \cos 2\phi +1 \right) d \phi = 0 + \frac{\pi}{2}.\]
Finally,
\[ T = 2T_e = \pi \sqrt{\frac{R_\text{max}^3}{2GM}} = \pi \sqrt{\frac{9R_\text{max}^3}{2c^2 R_\text{max}}} = \frac{3\pi R_\text{max}}{c\sqrt{2}}\]