You can observe that by Gauss Equation and Codazzi-Mainardi equation you can write the Riemann tensor and the Ricci tensor in function of the second fundamental form of your hypersurface that is in $\mathbb{R}^{n+1}$, so the Riemann and Ricci tensor of the ambient are zero.
The Weyl tensor is defined in function of the induced metric, the Riemann tensor and the Ricci tensor of the hypersurface so it is true, it is possible write the Weyl tensor in function of the second fundamental form.
We know that Gauss equation of an Hypersurface $M\subseteq N$ is
$R^N(X,Y,Z,W)=R^M(X,Y,Z,W)+h(X,W)h(Y,Z)-h(X,Z)h(Y,W)$
where $h$ is the second fundamental form of $M$.
In local coordinates we have that
$R^N_{ijkl}=R^M_{ijkl}+h_{il}h_{jk}-h_{ik}h_{jl}$
In our case $R^N=0$ so
$R^M_{ijkl}= h_{ik}h_{jl}-h_{il}h_{jk}$
If we trace the equation we can write also Ricci Tensor in function of the second fundamental form:
$R^M_{jl}=Hh_{jl}-g^{ik} h_{il}h_{jk}$
where $H$ is the mean curvature of $M$.
If we trace twice Gauss equation we get a relationship between the scalar curvature $R$ of $M$ and the second fundamental form:
$R=H^2-|h|^2$
We consider the (0,4)-Weyl tensor of $M$
$W=R^M-\frac{1}{n-2}(Ric-\frac{R}{n}g)\wedge g-\frac{R}{2n(n-1)}g\wedge g$
where $\wedge$ is the Kulkarni-Nomizou product.
In local coordinate we have that
$W_{ijkl}=R^M_{ijkl}+\frac{1}{n-2}(R^M_{il}g_{jk}-R^M_{ik}g_{il}+R^M_{jk}g_{il}-R^M_{jl}g_{ik})+$
$+\frac{R}{(n-1)(n-2)}(g_{ik}g_{jl}-g_{il}g_{jk})$
So you have all elements to write Weyl tensor in function of the second fundamental form of $M$ (and in function of the metric $g$).
You can observe also that the (0,2)-Tensor $h$ is symmetric so there exists a local orthonormal frame $\{Z_1,\dots, Z_n\}$ such that the matrix $A$ associated with $h$ with respect that frame is $diag(\lambda_1,\dots, \lambda_n)$ and th matrix $G$ associated with the metric $g$ with respect that frame will be the identity $I$. In this case you can write all our formula in an easier way than before.
If it is not clear I can explain more precisely.
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$.
Best Answer
You're correct that this is absurd; the error is in your definition of sectional curvature, discussed below. The computation about principal curvatures is correct, and gets at some interesting phenomena exhibited by minimal submanifolds. For instance, all embedded minimal $2$-manifolds of Riemannian manifolds are "saddle-shaped" in the sense that $\lambda_1 \lambda_2 \leq 0$, and as a consequence their sectional curvatures are bounded above by the sectional curvatures of the ambient space.
In the case of embedded $2$-dimensional submanifolds $(M^2, g)$ of $\mathbb{R}^3$, the scalar curvature of $M$ at $p$ is (twice) the product of the principal curvatures at $p$, which is independent of the realization of $(M^2,g)$ as a Riemannian submanifold of $\mathbb{R}^3$ by Gauss's Theorem Egregium. In this case, the observation is fully correct: every embedded minimal surface in $\mathbb{R}^3$ has nonpositive scalar curvature everywhere.
In general though, the products of principal curvatures of a Riemannian submanifold $(M,g) \subset (M',g')$ depend heavily on the realization of $(M,g)$ inside $(M',g')$, and upon $(M',g')$. Accordingly, the sectional curvatures of $M$ are defined purely in terms of $(M,g)$, and not in terms of an embedding of $(M,g)$ elsewhere.