Please consider the following metric
$$ ds^2 =-dt^2+f(t)( {\frac {a^2}{r^2+a^2}}dr^2+ r^2 d\Omega_{n-1})$$
where $d\Omega_{n-1}$ is the metric of the hyper-sphere $S^{n-1}$
Making "experimental symbolic computation" with the package GrTensor, I am obtaining that
$$f(t) = {\frac {1}{8}}{\frac {\left( 2\,{{\rm e}^{-2\,\sqrt {2}\sqrt {{\frac {\Lambda}{n \left( n-1 \right) }}}t}}-\frac{1}{2}\,n \left( n-1 \right) \right) ^{2}{{\rm e}^{2\,\sqrt {2}\sqrt {{\frac {\Lambda}{n \left( n-1 \right) }}} t}} }{{\Lambda}^{\frac{3}{2}}{a}^{3}}} $$
is a solution for the Einstein equation with cosmological constant:
$$ G_{{\mu,\nu}}+\Lambda\,g_{{\mu,\nu}}=0 $$
My question is: how to derive this solution using only pen and paper?
Possibly a full solution, but too long for a comment.
It would be much easier to solve this using the Friedmann equation in density parameter form. In $3+1$ dimensional space we have
$$\frac{\dot{f}^2}{H_0^2f^2} = \frac{1-\Omega_{0,\Lambda}}{f^2}+\Omega_{0,\Lambda}$$
since in your equation the universe is devoid of everything except universal constant and curvature. This equation is separable - and taking the solution branch where the universe is growing, we obtain
$$H_0(t-t_0)+C = \int\frac{df}{\sqrt{1-\Omega_{0,\Lambda}+\Omega_{0,\Lambda}f^2}} = \frac{1}{\sqrt{\Omega_{0,\Lambda}}}\sinh^{-1}\left(f\sqrt{\frac{\Omega_{0,\Lambda}}{1-\Omega_{0,\Lambda}}}\right)$$
$$\implies f(t) = \sqrt{\frac{1-\Omega_{0,\Lambda}}{\Omega_{0,\Lambda}}}\sinh\left(H_0\sqrt{\Omega_{0,\Lambda}}(t-t_0)+\sinh^{-1}\left(\sqrt{\frac{\Omega_{0,\Lambda}}{1-\Omega_{0,\Lambda}}}\right)\right)$$
or
$$f(t) = \cosh\Bigr(H_0\sqrt{\Omega_{0,\Lambda}}(t-t_0)\Bigr)+\frac{1}{\sqrt{\Omega_{0,\Lambda}}}\sinh\Bigr(H_0\sqrt{\Omega_{0,\Lambda}}(t-t_0)\Bigr)$$
This solves the relationship in $3+1$ dimensions but you have an $n+1$ dimensional space. This would require determining the Friedmann equations in $n+1$ dimensional space. I conjecture that the equation for a universal constant and curvature only universe would not actually change in this setting since the dimensionality would only affect the $T_{\mu\nu}$ side of the equation which is $0$ here. But I am not entirely sure.
If my conjecture is correct, then this is the full solution.