Calculating scalar curvature of a warped product $I \times N^n(k)$

My answer is below, but I wanted to make several comments.

First, regarding the ODE you found, the first two terms are $$ \frac{n}{u}\left(\frac{f'}{f} - ff'\right). $$ This is an immediate sign that you might have made an error, because it is unlikely you would get two terms that depend only $f$ and $f'$ but scale differently if you rescale $f$ by a constant factor. This is a hint to check particularly carefully the calculations that led to these two terms.

Second, I don't know if you tried rescaling the metric $g$. The observation is that, given $\alpha > 0$, if you set \begin{align*} f_\alpha &= \alpha^{-1}f\\ h_\alpha &= \alpha^2h\\ g_\alpha &= \frac{du^2}{f_\alpha^2} + u^2h_\alpha, \end{align*} then $g_\alpha = \alpha^2 g$. Therefore, if $S$ is the scalar curvature of $g$, then the scalar curvature of $g_\alpha$ is $$ S_\alpha = \alpha^{-2}S. $$ This means that if $f, h, S$ satisfy the ODE, then $f_\alpha, h_\alpha, S_\alpha$ had better satisfy the ODE, too. In your ODE, the second term scales properly, but the first and third do not. So that tells you not only that there is definitely an error but also where to focus your attention when searching for your error.

I know that error checking methods like this are rarely taught explicitly, especially not in textbooks, but developing such skills is quite important when doing differential geometric calculations.

The calculation below assumes you are familiar with using the moving frames method with differential forms. There are at least 3 ways to do calculations in differential geometry: 1) Using vector fields; 2) Using local coordinates; 3) Using moving frames. I recommend practicing and trying to develop facility with all 3. I find that depending on the calculation, often one way is easier than the other two. Even if you did not use moving frames, you can compare your formulas for the Riemann curvature, Ricci curvature, and scalar curvature against mine.

In general, I prefer to break down a calculation into natural pieces like this and do each calculation myself, instead of using a single formula that someone else has written down. There are a few reasons for this. That person might have made a mistake in their calculation or, even if they have the correct formula, copied it incorrectly. Or they are using different conventions than you are for the definitions and notation. Also, if I make a mistake using the formula, it is harder to identify where the error is than a calculation that is broken down into conceptually natural pieces.

Finally, here is my calculation, which I had to redo about 4 or 5 times. My final answer agrees with Robert's. That does not guarantee it is right, but it makes the probability high, especially since Robert is very reliable.

Let $\theta^1, \dots, \theta^n$ be an orthonormal frame of $1$-forms for the metric $h$ that has constant sectional curvature $\kappa$. Let $\theta^i_j$, $1 \le i,j \le n$, the corresponding connection $1$-forms. Therefore, \begin{align*} d\theta^i + \theta^i_j\wedge\theta^j &= 0\\ d\theta^i_j + \theta^i_p\wedge\theta^p_j &= \kappa\,\theta^i\wedge\theta^j. \end{align*} Given a function $f(u)$, consider the metric $$ g = \frac{du^2}{f^2} + u^2h. $$ An orthonormal frame of $1$-forms for $g$ is given by \begin{align*} \omega^0 &= \frac{du}{f}\\ \omega^i &= u\,\theta^i. \end{align*} Differentiating, \begin{align*} d\omega^0 &= 0\\ d\omega^i &= du\wedge\theta^i + u\,d\theta^i\\ &= -f\theta^i\wedge\omega^0 - \theta^i_j\wedge\omega^j \end{align*} This implies that the connection $1$-forms for $g$ with respect to the frame $(\omega^0, \dots, \omega^n)$ are \begin{align*} \omega^i_0 &= f\theta^i\\ \omega^i_j &= \theta^i_j. \end{align*} The curvature $2$-forms are therefore \begin{align*} \Omega^i_0 &= d\omega^i_0 + \omega^i_j\wedge\omega^j_0\\ &= f'\,du\wedge\theta^i + f\,d\theta^i + f\theta^i_j\wedge\theta^j\\ &= \frac{-ff'}{u}\,\omega^i\wedge\omega^0\\ \Omega^i_j &= d\omega^i_j + \omega^i_0\wedge\omega^0_j + \omega^i_p\wedge\omega^p_j\\ &= d\theta^i_j + \theta^i_p\wedge\theta^p_j - f^2\theta^i\wedge\theta^j\\ &= \frac{\kappa^2-f^2}{u^2}\omega^i\wedge\omega^j. \end{align*} This implies that \begin{align*} R_{i0i0} &= \frac{-ff'}{u}\\ R_{ijij} &= \frac{\kappa-f^2}{u^2}, \end{align*} where $1 \le i,j\le n$ and $i \ne j$, and all other components are zero. It follows that the only possible nonzero components of the Ricci curvature tensor are \begin{align*} R_{00} &= \sum_{i=1}^n R_{i0i0} = \frac{-nff'}{u}\\ R_{ii} &= R_{i0i0} + \sum_{j\ne i} R_{ijij}\\ &= -\frac{ff'}{u} + \frac{(n-1)(\kappa-f^2)}{u^2}. \end{align*} This implies that the scalar curvature is \begin{align*} S &= R_{00} + \sum_{i=1}^n R_{ii}\\ &= \frac{-nff'}{u} + n\left(-\frac{ff'}{u} + \frac{(n-1)(\kappa-f^2)}{u^2}\right)\\ &= \frac{-2nff'}{u} + \frac{n(n-1)(\kappa-f^2)}{u^2}\\ &= \frac{n}{u^2}(-2ff'u + (n-1)(\kappa-f^2)). \end{align*} As written, it is not clear how to solve this ODE. However, it is easy to see that since $2ff' = (f^2)'$, you can rewrite the equation in terms of $$\ \phi = \kappa-f^2. $$ The equation can be rewritten as $$ u\phi' + (n-1)\phi = \frac{S}{n}u^2 $$

If we now multiply by $u^{n-2}$, then we get $$ u^{n-1}\phi' + (n-1)u^{n-2}\phi = \frac{S}{n}u^{n}, $$ which implies $$ (u^{n-1}\phi)' = \frac{S}{n}u^{n} $$ If we assume that $S$ is constant and $$ \lim_{u\rightarrow 0} u^{n-1}\phi(u) = c, $$ then the fundamental theorem of calculus implies that $$ u^{n-1}\phi = \frac{S}{n(n+1)}u^{n+1} + c. $$ Solving for $\phi$, we get $$ \phi = \frac{S}{n(n-1)}u^2 + cu^{1-n}. $$ Therefore, \begin{align*} f^2 &= \kappa - \phi\\ &= \kappa - \frac{S}{n(n+1)}u^2 - cu^{1-n}. \end{align*} It follows that the metric $$ g = \frac{du^2}{\kappa- \frac{S}{n(n+1)}u^2 - cu^{1-n}} $$ has constant scalar curvature $S$.