Is it possible to compute explicitly the integral curves of the vector field: $$X=(y^2+z^2)\,\partial_x-xy\,\partial_y+xz\,\partial_z$$ The first integrals of $X$ are $yz$ and $x^2+y^2-z^2$ but its integral curves are hard to compute.
Do you have any suggestion?
Here is a suggestion (but not a complete solution). The integral curves are the solutions to the system of differential equations \begin{align} \dot{x}(t)&=y^2+z^2, \tag{1a} \\ \dot{y}(t)&=-xy, \tag{1b} \\ \dot{z}(t)&=xz, \tag{1c} \end{align} with initial conditions $x(0)=x_0, y(0)=y_0, z(0)=z_0$. Now, as you already noted, $x^2+y^2-z^2$ is a first integral of $(1)$, i.e., $$ x(t)^2+y(t)^2-z(t)^2=C, \tag{2} $$ where $C(=x_0^2+y_0^2-z_0^2)$ is a constant.
Using $(2)$, it's possible to decouple the system $(1)$. Indeed, \begin{align} \ddot{x}&=2y\dot{y}+2z\dot{z}=-2yxy+2zxz=2x(-y^2+z^2) =2x(x^2-C), \tag{3a} \\ \ddot{y}&=-\dot{x}y-x\dot{y}=-(y^2+z^2)y+x^2y=(C-2y^2)y, \tag{3b} \\ \ddot{z}&=\dot{x}z+x\dot{z}=(y^2+z^2)z+x^2z=(C+2z^2)z. \tag{3c} \end{align} Since the equations $(2)$ are of second order, they need two initial conditions each. The extra initial conditions can be derived from $(1)$: $\dot{x}(0)=y_0^2+z_0^2, \dot{y}(0)=-x_0y_0, \dot{z}(0)=x_0z_0.$
Equations $(3)$ can be reduced to quadratures. For instance, multiplying $(3\text{a})$ by $\dot{x}$ and integrating, one obtains $$ \frac{1}{2}\dot{x}^2=\frac{1}{2}x^4-Cx^2+C_1 \implies t=\int_{x_0}^x\frac{du}{\sqrt{u^4-2Cu^2+2C_1}}, \tag{4} $$ where $C$ and $C_1$ are functions of the initial conditions, and the sign of the square root is positive because, according to $(1\text{a})$, $\dot{x}(t)\geq 0$. Similar integrals can be derived from $(3\text{b})$ and $(3\text{c})$.
Addendum
Simpler expressions can be achieved by reparametrizing time. Let's define a new time variable $\tau$ satisfying the differential equation $$ \dot{\tau}(t)=xz, \qquad \tau(0)=0. \tag{A.1} $$ Combining $(\text{A}.1)$ with $(1\text{c})$ we obtain $$ \frac{dz}{d\tau}=\frac{\dot{z}}{\dot{\tau}}=1 \implies z(\tau)=\tau+z_0. \tag{A.2} $$ From $(\text{A}.1)$ and $(1\text{b})$ it follows that $$ \frac{dy}{d\tau}=\frac{\dot{y}}{\dot{\tau}}=-\frac{y}{\tau+z_0} \implies y(\tau)=\frac{y_0z_0}{\tau+z_0}. \tag{A.3} $$ (This result can also be obtained from $(\text{A}.2)$ and the first integral $yz=y_0z_0$.)
Finally, combining all the above results with $(1\text{a})$, we find $$ \frac{dx}{d\tau}=\frac{\dot{x}}{\dot{\tau}}=\frac{y^2+z^2}{xz} \implies x\frac{dx}{d\tau}=\frac{y_0^2z_0^2}{(\tau+z_0)^3}+\tau+z_0 $$ $$ \implies x(\tau)^2=x_0^2+y_0^2-\frac{y_0^2z_0^2}{(\tau+z_0)^2}+\tau^2+2z_0\tau. \tag{A.4} $$ (This result can also be obtained directly from $(2)$, $(\text{A}.2)$, and $(\text{A.3})$.)