Existence of a solution for a quasilinear hyperbolic system of PDEs with many state variables I'll star by saying that I am not really familiar with the field of PDEs so this questions may be trivial or ill-possed in that case please let me know.
I am  in search of some existence (Global) result regarding a system of first order PDEs with many state-variables and a non-homogeneity that is non-linear in the solutions, i.e.
\begin{align}
\begin{cases}
\partial_t \mathbf u+\sum_{n=1}^{K} A_n(t)\partial_{x_n} \mathbf u=\mathbf b(t,\mathbf u,x_1,\cdots,x_K)\\
\mathbf u(0,x_1,\cdots,x_K)=u_0\in\mathbb R^N,
\end{cases}
\end{align}
where the matrices $A_n, n=1,\cdots, K$ are diagonal, and $\mathbf b$ is Lipschitz continuous in $\mathbf u$ uniformly in $t,x_1,\cdots,x_n$.
In Bressan's "Conservations laws" book the case with two variables is considered, namely $t$ and $x$. I have found some articles online considering the case with many state variables but in that case the non-homogeneity was linear in the $\mathbf u$ .
As far as I know the method of the characteristics applies only the one-dimensional case, so I am a little bit lost.
Is there some result I can use to show the existence (and uniqueness if possible) of a solution for this system?
Thanks in advance.
 A: Okay, so I would write your equations instead in the following form:
$$ \partial_t u_i + v_i(t) \cdot \nabla u_i = b_i (t, \vec{x}, \vec{u}) $$
This is a system of transport equations and so can actually be solved by using a variation of the Picard-Lindelof argument.
(I am implicitly assuming that your function $b$ is suitably nice in certain ways, which are implied by what you said is allowed in this comment.)
1
Consider the iteration scheme
$$ \vec{w} \mapsto \vec{w}' $$
where $\vec{v}'$ solves the linear inhomogeneous equation
$$ \partial_t w'_i + v_i(t) \cdot \nabla w'_i = b_i(t,\vec{x}, \vec{w})$$
with the fixed, given, initial data $\mathring{\vec{w}}:\mathbb{R}^K\to \mathbb{R}^N$.
Given a continuous function $\vec{w}: [-t_0, t_0] \times\mathbb{R}^K \to \mathbb{R}^N$, the function $\vec{w}'$ can be solved component wise by integrating along the integral curves of $\partial_t + v_i(t) \cdot \nabla$ on $[-t_0, t_0]\times \mathbb{R}^K$.
It thus suffices to show that this iteration scheme converges.
2
Suppose $\vec{z}$ and $\vec{w}$ are given functions that are uniformly bounded by some constant $M$. Then there exists (under what I hope are reasonable assumptions on $b_i$) some constant $M'$ such that $b_i(t,\vec{x},\vec{z})$ and $b_i(t,\vec{x},\vec{w})$ are uniformly bounded by $M'$. Thus you have that
$$ | \vec{z}'|_\infty, |\vec{w}'|_\infty \leq |\mathring{\vec{w}}|_\infty + M' t_0 $$
Choose $M = 2 |\mathring{\vec{w}}|_\infty$ and there exists some $t_0 > 0$ such that the iteration mapping maps the (closed) ball of radius $M$ in $C^0([-t_0, t_0]\times \mathbb{R}^K,\mathbb{R}^N)$ to itself.
3
For the differences you find that
$$ \partial_t (z'_i - w'_i) + v_i(t)\cdot \nabla (z'_i - w'_i) = b_i(t,\vec{x}, \vec{z}) - b_i(t,\vec{x},\vec{w}) $$
Using the Lipschitz continuity in the final slot of $b_i$ you get that
$$ \big|\partial_t (z'_i - w'_i) + v_i(t)\cdot \nabla (z'_i - w'_i) \big| \leq M'' |\vec{z} - \vec{w}| $$
for some $M''$. Integrating you get that if $t_0$ is chosen sufficiently small, the iteration mapping is a contraction mapping.
4
Applying Banach fixed-point theorem you get local existence.
5
Having proven local existence, you can upgrade this to global in the same way you argue for ODEs: by showing that the sup norm of the functions involved does not blow-up. This follows from the global Lipschitz property of the functions $b_i$, which guarantees that $|b_i(t,\vec{x},\vec{w})| \leq M''' (1 + |\vec{w}|)$
So a standard Gronwall's inequality argument well tell you that the sup-norm of the solution cannot grow faster than $A e^{2 M''' |t|}$ for some constant $A$.
