You get a separate ODE for every integral curve (which you denote by $\gamma$) of the vector field $X$. The space $\mathbb R^4$ is a disjoint union of these integral curves. Each ODE is of order one, so you need to fix the value of $f$ at exactly one point from every integral curve. (Alternatively, you could also take asymptotic behaviour on each integral curve as a boundary condition.) You can fix the values on a three dimensional surface that meets every integral curve exactly once, provided such a surface exists. A good choice for such a surface or a collection of such surfaces depends on the vector field $X$.
The curves $\gamma$ should be curves in $\mathbb R^4$. We want $\gamma$ to satisfy $\gamma'(t)=X(\gamma(t))$ for every $t\in\mathbb R$ (this is what an integral curve mean) and we want $f$ to satisfy the ODE $\frac{d}{dt}f(\gamma(t))=u(\gamma(t))f(\gamma(t))$ along $\gamma$ (this for all integral curves $\gamma$ is equivalent with $Xf=uf$ if $f\in C^1$). This ODE can be solved by standard techniques.
The exact choice of the hypersurface on which we want to fix the values of $f$ depends on the vector field $X$. If $X(x)=(0,1+x_3^2,0,0)$, then one possibility is to fix the values of $f$ on the hypersurface given by $x_2=0$. There is no choice of the hypersurface that works for all possible $X$.