We consider the first-order Cauchy problem \documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$
\begin{gathered}
\partial _z u + a(z,x,D_x )u = 0,0 < z \leqslant Z, \hfill \\
u|_{z = 0} = u_0 , \hfill \\
\end{gathered}
$$\end{document} with Z > 0 and a(z, x,Dx) a k × k matrix of pseudodifferential operators of order one, whose principal part a1 is assumed symmetrizable: there exists L(z, x, ξ) of order 0, invertible, such that \documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$
a_1 (z,x,\xi ) = L(z,x,\xi )( - i\beta _1 (z,x,\xi ) + \gamma _1 (z,x,\xi ))(L(z,x,\xi ))^{ - 1} ,
$$\end{document} where β1 and γ1 are hermitian symmetric and γ1 ≥ 0. An approximation Ansatz for the operator solution, U(z′, z), is constructed as the composition of global Fourier integral operators with complex matrix phases. In the symmetric case, an estimate of the Sobolev operator norm in L((H(s)(Rn))k, (H(s)(Rn))k) of these operators is provided, which yields a convergence result for the Ansatz to U(z′, z) in some Sobolev space as the number of operators in the composition goes to ∞, in both the symmetric and symmetrizable cases. We thus obtain a representation of the solution operator U(z′, z) as an infinite product of Fourier integral operators with matrix phases.