Solución parcial. Es fácil ver que $ f $ tiene una inversa $ g $. Vamos a resolver el problema con la hipótesis de que la $ g $ es una función polinómica.
Por el teorema de la función inversa, ya que $f^\prime(x)> 0$ para todos los $ x \in [0, \infty) $ entonces $f:[0,\infty)\to [f(0),\infty)$ tiene una inversa global $g:[f(0),\infty)\to [0,\infty)$. Aviso de que se considere la posibilidad de $ f^\prime(0) $ como es derivable en a la derecha de cero. Además tenemos las igualdades
$$
\begin{matrix}
g(f(x))=x \quad & \quad f(g(y))=y \quad & \quad g(y)=x \quad & \quad f(x)=y\\
\\
g^\prime(y)=\frac{1}{f^\prime (x)} \quad & \quad g^\prime(y)=\frac{1}{f^\prime (g(y))}
\quad & \quad g(s)=t \quad & \quad f(t)=s
\end{de la matriz}
$$
Por la fórmula de cambio de variables que hemos
\begin{align}
\int_0^{t} \big( f(x) \big)^{1+\epsilon}\frac{1}{f^\prime(x)} \mathrm{d}x
=&
\int_{g(f(0))}^{g(s)} \big( f(x) \big)^{1+\epsilon}\frac{1}{f^\prime(x)} \mathrm{d}x,
\\
=&
\int_{f(0)}^{s} \big( f(g(y)) \big)^{1+\epsilon}\frac{1}{f^\prime(g(y))}g^\prime(y) \mathrm{d}y,
\\
=&
\int_{f(0)}^{s} \big( y \big)^{1+\epsilon}\cdot g^{\prime}(y)\cdot g^{\prime}(y)\mathrm{d}y,
\end{align}
Ahora, de verificación para cada función polinómica $ g(y)=a_ny^n+\ldots+a_1 y+a_0$que
$$
\Phi_g(s)=
\frac{1}{g(s)^2}\int_{f(0)}^{s} \big( y \big)^{1+\epsilon}\cdot g^{\prime}(y)\cdot g^{\prime}(y)\mathrm{d}y,
=
\frac{1}{t^2}\int_0^{t} \big( f(x) \big)^{1+\epsilon}\frac{1}{f^\prime(x)} \mathrm{d}x
$$
es una función tal que $\lim_{s\to\infty}\Phi_{g}(s)=\infty$.