There may not be a solution defined in $[0,\infty)$. Consider for each $a\in\mathbb{R}^+$ the curve $\gamma_a:[0,\frac{\pi a}{2})\to\mathbb{R}^2$ given by $\gamma_a(t)=(t,\text{tan}(\frac{t}{a})-ta)$, for $t\in[0,\frac{\pi a}{2})$. E.g., here are the curves $\gamma_a$ for $a=0.6,0.9,1.2,1.5$, in red, green, blue and orange respectively.
We will show that there is an initial value problem of the requested form for which the $\gamma_a$ are the only solutions, thus providing solutions on any $[0,b)$ but not on $[0,\infty)$.
The differential equation will use an appropriate inverse to the above function $\gamma$, so we begin by defining that inverse.
Each point $(t,x)\in(0,\infty)\times\mathbb{R}$ belongs to $\gamma_a$ for exactly one value of $a$, because such $a$ has to satisfy $t<\frac{\pi a}{2}$, so $a>\frac{2t}{\pi}$, and $x=\tan(\frac{t}{a})-ta$, and the function $a\mapsto\tan(\frac{t}{a})-ta$ decreases from $\infty$ to $-\infty$ when $a$ goes from $\frac{2t}{\pi}$ to $\infty$. Also, $a$ depends continuously on the point $(t,x)$.
In fact, $a$ depends smoothly on $(t,x)$: this is because, near a given point $(x_0,t_0)$, the values of $a$ are given by the function $\phi(x,t)$ such that, if $F(x,t,a)=x-\tan(\frac{t}{a})+ta$, then we have $F(x,t,\phi(x,t))=0$. And $\frac{\partial F}{\partial a}=t+\frac{t}{\cos^2(\frac{t}{a})a^2}\neq0$, so by the implicit function theorem $\phi$ is some $C^\infty$ function of $x,t$.
So if we arrange the function $f:(0,\infty)\times\mathbb{R}\to\mathbb{R}$ so that the solutions of the differential equation follow the curves $\gamma_a$, then none of these solutions are defined on $[0,\infty)$. Such a function $f$ can be expressed explicitly as $f(t,x(t))=-a-\frac{1}{\cos^2(\frac{t}{a})a}$, where $a=\phi(x,t)$. Indeed, this guarantees that for any fixed $a$, the curves $x(t)=\tan(\frac{t}{a})-ta$ are solutions of $x'(t)=f(x,t)$ for $t\in(0,\frac{\pi a}{2})$.
Moreover, those are all the solutions to the equation. This is because if a solution $x(t)$ passes through some point $(x_0,t_0)$ with $t_0>0$, then the curve $\gamma(t)=(t,x(t))$ satisfies the equation $\gamma'(t)=(1,f(x,t))$. So as the function $(x,t)\mapsto(1,f(x,t))$ is $C^\infty$, by the Picard Lindelöf theorem we have a unique solution $\gamma(t)=\gamma_a(t)$ for all $t\in(0,\frac{\pi a}{2})$, where $a=\phi(t_0,x_0)$.
No. Consider $x'=g(x,y)$, $y'=h(x,y)$. If we take $g(x,y)=2|x|^{1/2}$ for $y=x^2$ similarly $h=4|x|^{3/2}$ on $y=x^2$, then we
can check directly that $x=t^2$, $y=t^4$, $t\ge 0$, and $x=y=0$ are solutions with initial value $(0,0)$, independently of how we define $g,h$ off the parabola.
To define $g,h$ everywhere, we interpolate linearly, starting from the value on the parabola, by letting $g(ta,ta^2)=tg(a,a^2)$, $a\in\mathbb R$, $-1\le t\le 1$, and we also set $g(x,0)=g(0,y)=0$.
So far, $g$ has been defined on the closed set $\{ |y|\ge x^2\}\cup \{ (x,0)\}$, and it is continuous there.
To verify continuity on the $y$ axis, note that if $(ta,ta^2)\to (0,b)$ with $b\not= 0$, then $|a|\to\infty$, so $tg(a,a^2)=2t|a|^{1/2}\to 0$ (and also $t|a|^{3/2}\to 0$, which we need in the corresponding step for $h$).
The existence of all directional derivatives at the origin is already guaranteed because the current domain contains a line segment in every direction, and $g$ is linear there. Finally, we extend $g$ to $\{ 0<|y|<x^2\}$. We can make $g$ smooth away from the origin.
Best Answer
This is a partial solution, assuming that $f$ is analytic on $I\times R$, and so solutions are analytic as well, and that the space $E$ of solutions is of the same dimension $n$ as vectors $y$ and $f$.
We represent solutions $y$ and the function $f$ in the right hand side as column vectors, as usual. Let $y_1,\ldots,y_n$ be a basis in $E$. Then every vector $y\in E$ can be written as a linear combination $$y=c_1y_1+\ldots+c_ny_n.\quad\quad\quad\quad\quad\quad\quad\quad (1)$$ Solving this by Cramer's rule we obtain $c_k$ as a ratio of two determinants made of coordinates of $y$ and $y_j$. Since the coordinates of $y$ make one column in the numerator, it is clear that this ratio is a linear function of coordinates of $y$, so we can write $$c_j=a_j(t)y,\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad (2)$$ for some row vector $a_j$ which depends on $t$. Now, since all elements of $E$ satisfy your non-linear equation we must have for every $y\in E$ of the form (1): $$y'=\sum_{j=1}^n c_jy_j^\prime=\sum_{j=1}^nc_jf(t,y_j),$$ Substituting our formula for $c_j$ we obtain that $$y'=A(t)y,\quad\mbox{where}\quad A(t)=\sum_jf(t,y_j(t))a_j(t)$$ (column vectors $f(t,y_j(t))$ times row vectors $a_j$ give you $n\times n$ matrices.) There is a little problem with the denominator in the Cramer Rule: one has to show that it cannot be zero at some point. If all functions $f,y_j$ are analytic, this can only happen at isolated points. Away from these isolated points $t_k$ both equations $y'=f(t,y)$ and $y'=A(t)y$ have the same set of solutions, therefore $f(t,y)=A(t)y$ on some open set of $(t,y)$, and since all functions are analytic they must coincide.