What lhf says in the comments is generally true: Most of the time, there is no hope of finding a closed form solution. However, in this case, we can do it.
Let $\gamma(t) = (\gamma_1(t), \gamma_2(t), \gamma_3(t)$ be an integral curve with initial point $\gamma(0) = (x_0, y_0, z_0)$. What does this mean?
It means that $u$ at the point $\gamma(t)$ is equal to \gamma'(t). Let's write this out. I'm going to use $\partial_k$ for $\dfrac{\partial}{\partial x_k}$ to save typing.
$u(\gamma(t)) = \gamma_1 \partial_2 + \gamma_2 \partial_1 + \gamma_3\partial_3$ and \gamma'(t) = (\gamma_1'(t), \gamma_2'(t), \gamma_3'(t)) = \gamma_1' \partial_1 + \gamma_2'\partial_2+\gamma_3' \partial_3.
Setting these equal to each other and equating coefficients gives us a system of ODEs to solve:
\gamma_1' = \gamma_2 \gamma_2' = \gamma_1 \gamma_3' = \gamma_3
The first two equations are coupled but the third is not, so lets solve that one first. The solution to \gamma_3' = \gamma is $\gamma(t) = Ce^t$ for some constant $C$.
There is a known process for solving coupled linear ODEs, but in this case, I think it's easier to just guess a solution. We want two functions so that if we start with one and take two derivatives, we get back where we started. This suggests we try $\gamma_1(t) = Ae^t + Be^{-t}$. Plugging this into the second equation gives $\gamma_2(t) = Ae^t-Be^{-t}$, and it's easy to check that this choice of $\gamma_1$ and $\gamma_2$ solves the second equation.
The upshot is we now know $\gamma(t) = (Ae^t + Be^{-t}, Ae^t -Be^{-t}, Ce^t)$. What are $A$, $B$, and $C$?
Well, $\gamma(0) = (A+B,A-B, C) = (x_0,y_0,z_0)$. So, $A = \frac{x_0+y_0}{2}$ and $B = \frac{x_0-y_0}{2}$ while $C = z_0$.