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]
