Cfrgtkky

Consider the following system of partial differential equations defined on some neighborhood $U$ of the origin in $mathbb{R}^2:$
$$z_x = F(x,y,z(x,y)), z_y = G(x,y,z(x,y))$$ where $F$ and $G$ are smooth on $U times mathbb{R}.$

In an exercise in the book "Manifolds and Differential Geometry" by Jeffrey M. Lee it is claimed that solving the above system of differential equations is equivalent to finding the integral curves of the pair of vector fields $$X = dfrac{partial}{partial x} + F(x,y,z)dfrac{partial}{partial z},$$ $$Y = dfrac{partial}{partial y} + G(x,y,z) dfrac{partial}{partial z}.$$

I do not see why this is true. Clearly, finding an integral curve for $X$ is equivalent to finding a curve $gamma(t)= (gamma_1(t),gamma_2(t),gamma_3(t))$ such that $$gamma_1'(t) = 1,$$ $$gamma_2'(t) = 0,$$ $$gamma_3'(t) = F(gamma_1(t),gamma_2(t),gamma_3(t)).$$ One could imagine by putting some certain initial conditions that solutions are given by $gamma_1(t) = x,$ $gamma_2(t) = a$ for $a$ some constant, so that the last equation reads
$$gamma_3'(t) = F(t,a,gamma_3(t)).$$ Finding a solution to this equation should then be thought of as solving the differential equation $z_x(t,a) = F(t,a,z(t,a)).$ We can clearly put these solutions together using a flow $varphi^X_(t,a).$ An analogous analysis can then be made for finding integral curves to $Y.$

Here is my question:
How do I combine these individual solutions to a global solution?
In fact, I do not see how to do this unless we have some more conditions on the vector fields, such as that $[X,Y]$ vanishes.