System of linear pde with non constant coefficients

If you switch the second and third rows of your system, the differential operator is the same as the linearization of equation (4.3) in Existence of elastic deformations with prescribed principal strains and triply orthogonal systems at a point where $$ g_{\alpha\beta} = \delta_{\alpha\beta} \text{ and }\frac{\partial x^\alpha}{\partial y^i} = \delta_{\alpha i}. $$ As mentioned immediately after the equation, this system is not (symmetric) hyperbolic. It also can never be elliptic. It follows that this system is difficult to solve. I don't know of any work in this direction.

In the paper cited, we observe that the obstruction is due to a type of gauge invariance and reformulate the problem in a way that eliminates this issue, resulting in a symmetric hyperbolic system that can be solved.

Perhaps your problem can be treated similarly?

Recall that a first order system $A_{ij}^\mu \partial_\mu \phi^j + B_{ij} \phi^j = 0$ (the right hand-side could also be inhomogeneous) is symmetric hyperbolic when there exists at least one covector $p_\mu$ such that $A_{ij}^\mu p_\mu > 0$, the contraction is a positive definite symmetric matrix. If the coefficients $A_{ij}^\mu = A_{ij}^\mu(x)$ are $x$-dependent, then the positivity condition should be satisfied for all $x$. There are no special conditions imposed on $B_{ij}$, which could also be $x$-dependent. This is all standard. And the initial value problem for such systems is well posed, at least on those domains where the causal structure of formed by the cones of $p_\mu$'s satisfying the above positivity conditions is globally hyperbolic (the domain admits a foliation by Cauchy surfaces). See for instance the excellent account in Chapter 7 of [1].

Let me now introduce some non-standard terminology. I call a linear system $$A_{ij}^{\mu_1\cdots \mu_k} \partial_{\mu_1} \cdots \partial_{\mu_k} \phi^j + B(\phi,\partial\phi, \cdots \partial^{k-1}\phi) = 0 \tag{*}$$ higher order symmetric hyperbolic when it can be put in symmetric hyperbolic form by reduction to first order. Next, I call a linear system (using the notation $A_{ij}(\partial) = A_{ij}^{\mu_1\cdots \mu_k} \partial_{\mu_1} \cdots \partial_{\mu_k}$ and writing l.o.t for what I wrote as $B_{ij}(\cdots)$ above) $$A_{ij}(\partial) \phi^j + l.o.t = 0$$ generalized symmetric hyperbolic when there exists a complementary operator $C_k^k(\partial) = (C_k^i)^{\mu_1\cdots \mu_l} \partial_{\mu_1} \cdots \partial_{\mu_l}$ such that applying it to (*) gives a system $$ C_k^i(\partial) A_{ij}(\partial) \phi^j + l.o.t = 0$$ that is higher order symmetric hyperbolic.

The point is that a generalized symmetric hyperbolic system has an equally well-posed initial value problem, with or without a source term, as a symmetric hyperbolic one. Unfortunately, I don't have a reference for this, except for the parallel discussion that I gave for generalized normally hyperbolic systems in the recent paper [2] (just replace "normally hyperbolic" by "symmetric hyperbolic" everywhere).

Claim: Your PDE system is generalized symmetric hyperbolic.

Let me first illustrate by a relevant example, what it means for a system to be higher order symmetric hyperbolic. Take the equation $$\delta_{ij} \partial_x \partial_y \partial_z N^j + l.o.t = 0 , \tag{**}$$ where $\delta_{ij}$ is just the Kronecker delta. By introducing the auxiliary variables $N_z^j = \partial_z N^j$ and $N_{yz}^j = \partial_y N_z^j$, it can be reduced to the first order system $$\begin{pmatrix} \delta_{ij}\partial_z & 0 & 0 \\ 0 & \delta_{ij}\partial_y & 0 \\ 0 & 0 & \delta_{ij}\partial_x \end{pmatrix} \begin{pmatrix} N^j \\ N_z^j \\ N_{yz}^j \end{pmatrix} + l.o.t = 0 ,$$ where the desired positivity is satisfied by any covector from the first octant, $p_x,p_y,p_z>0$. This shows that (**) is indeed higher order symmetric hyperbolic.

Your PDE system has the special form $$\begin{pmatrix} f_2 \partial_y & f_1 \partial_x & 0 \\ 0 & f_3 \partial_z & f_2 \partial_y \\ f_3 \partial_z & 0 & f_1 \partial_x \end{pmatrix} \begin{pmatrix} N_1 \\ N_2 \\ N_3 \end{pmatrix} = \begin{pmatrix} \omega_1 \\ \omega_2 \\ \omega_3 \end{pmatrix} .$$ As mentioned in Deane Yang's answer, this equation cannot be put into symmetric hyperbolic form. That is, there exists no matrix $C_k^i$ such that the principal symbol $C_k^i A_{ij}^\mu p_\mu$ satisfies both the symmetry and positivity conditions. You can find $C_k^i$ such that the symbol becomes symmetric, but not positive. At least one eigenvalue of the resulting $p$-dependent matrix will always remain negative.

However, assuming that $f_1,f_2,f_3 \ne 0$ are nowhere vanishing, one can multiply this system by a second order matrix differential operator (which you can read off from the formula below) such that the system becomes $$\partial_x \partial_y \partial_z \begin{pmatrix} N_1 \\ N_2 \\ N_3 \end{pmatrix} + l.o.t = \frac{1}{2 f_1 f_2 f_3} \begin{pmatrix} f_1 f_3 \partial_x \partial_z & -f_1^2 \partial_x^2 & f_1 f_2 \partial_x \partial_y \\ f_2 f_3 \partial_y \partial_z & f_1 f_2 \partial_x \partial_y & -f_2^2 \partial_y^2 \\ -f_3^2 \partial_z^2 & f_1 f_3 \partial_x \partial_z & f_2 f_3 \partial_y \partial_z \end{pmatrix} \begin{pmatrix} \omega_1 \\ \omega_2 \\ \omega_3 \end{pmatrix} .$$ From the example (**) it should now be obvious that this system is higher order symmetric hyperbolic, verifying the claim that your system is generalized symmetric hyperbolic. Therefore, it has a well-posed initial-value problem on any hyperplane whose conormal $p_\mu$ lies in the first octant.

[1] Ringström, Hans, The Cauchy problem in general relativity., ESI Lectures in Mathematics and Physics. Zürich: European Mathematical Society (EMS) (ISBN 978-3-03719-053-1/pbk). xiii, 294 p. (2009). ZBL1169.83003.

[2] See the discussion around Lemma 3 in García-Parrado, Alfonso and Khavkine, Igor, Conformal Killing Initial Data, Journal of Mathematical Physics 60 122502 (2019). [arXiv:1905.01231]