Non-homogeneous Linear System
I'm gonna consider $u = u(t)$ (and $b \in \mathbb{R}$) as a function of time. We let $\dot{x}(t) := \mathrm dx/\mathrm dt$ then
$$\left\{\begin{matrix}\dot{x}_1 = a_{11}x_1\\\dot{x}_2 = a_{21}x_1+a_{22}x_2 + bu(t)\end{matrix}\right.$$
Which can be transformed in a linear vector differential equation
$$\dot{X} = AX + U(t) \tag{1}$$
Where we have that
$$X := \left(\begin{matrix}x_1\\x_2\end{matrix}\right) \hspace{3cm}A := \left(\begin{matrix}a_{11}&0\\a_{21} & a_{22}\end{matrix}\right) $$
$$U(t):=\left(\begin{matrix}0\\bu(t)\end{matrix}\right)$$
We consider $A$ constant and $U = U(t)$ time-dependent. So we have to solve the equation $(1)$.
Existence and Uniqueness
Consider that a initial condition $X_0 := X(0)$ is given, then the equation $(1)$ has one and only one solution $X(t)$ that satisfies the initial condition. This is known as the Cauchy problem. The proof for that is so beautiful, in my opinion, that, even a little long, deserves to be shown. Let us first consider the homogeneous differential equation with the below Cauchy problem:
$$\dot{X} = AX \hspace{3cm}X(0):=0$$
Where $A$ is a $n\times n$ matrix and $X$ is a $n$-dimensional vector. If we integrate both sides of this equation we have that
$$X(t) = \int_0^tAX(t_1)\mathrm dt_1$$
Using this fact interactively for $m$ times we have we have
$$X(t) = \int_0^t\int_0^{t_1}\dots\int_0^{t_{m-1}}A^mX(t_m)\mathrm dt_m\dots\mathrm dt_1$$
If we look for the $i$-coordinate of the vector $X$ this equation is saying
$$x_i(t) = \int_0^t\int_0^{t_1}\dots\int_0^{t_{m-1}}\sum_{j=1}^n\left(\underbrace{A\cdot A\cdot A\dots A}_{\text{m times}}\right)_{ij}x_j(t_m)\mathrm dt_m \dots \mathrm dt_1$$
Where I used that
$$(AX)_i = \sum_{j=1}^nA_{ij}X_j = \sum_{j=1}^nA_{ij}x_j$$
For product of matrix we have that
$$(AB)_{ij} = \sum_kA_{ik}B_{kj}$$
Then a product of $m$ matrices we have $m-1$ sums so
$$\left(A^m\right )_{ij} = \sum_{k_1 = 1}^n\sum_{k_2 = 1}^n \dots \sum_{k_{m-1}=1}^nA_{ik_1}A_{k_1k_2}\dots A_{k_{m-1}j}$$
This implies that
$$x_i(t) = \sum_{j=1}^n\sum_{k_1 = 1}^n\sum_{k_2 = 1}^n \dots \sum_{k_{m-1}=1}^n\int_0^t\int_0^{t_1}\dots\int_0^{t_{m-1}}A_{ik_1}A_{k_1k_2}\dots A_{k_{m-1}j}x_j(t)\mathrm dt_m \dots \mathrm dt_1$$
This implies that the absolute value of $x_i(t)$ respects the inequality
$$\vert x_i(t)\vert \leq \sum_{j=1}^n\sum_{k_1 = 1}^n\sum_{k_2 = 1}^n \dots \sum_{k_{m-1}=1}^n\int_0^t\int_0^{t_1}\dots\int_0^{t_{m-1}}\vert A_{ik_1}\vert \vert A_{k_1k_2}\vert \dots\vert A_{k_{m-1}j}\vert \vert x_j(t)\vert \mathrm dt_m \dots \mathrm dt_1$$
Now suppose (in this moment) that $t \in [0,T]$ for some finite $T \in \mathbb{R}$ then we define $\alpha := \max_{i,j\in\{1,\dots,n\}}\vert A_{ij} \vert$ and $M := \max_{t \in [0,T]}\max_{i \in \{1,\dots,n\}}\vert x_i(t)\vert$ then the inequality becomes
$$\vert x_i(t) \vert \leq \alpha^mM\frac{t^m}{m!}n^m$$
Then we have that this inequality must hold for $t \in [0,T]$ for every $i$-component of $X$. Note that the right-side of the equation doesn't depends on $m$, this is just the number of interaction made. If we set $m \rightarrow \infty$ we have that $x_i(t)=0$ for $t \in [0,T]$ for all $i$ and for all $t$. We set $T$ arbitrarily so this is true for all $t$ and $i$. We conclude that the only solution is the vector $X = 0$ for the homogeneous equation with initial condition $X_0 = 0$. Now take equation $(1)$ and let $X_1$ and $X_2$ be two different solutions of equation $(1)$ with initial condition $X_1(0) = X_2(0) = X_0$. So let define $Z := X_1 - X_2$ so we have that $Z(0)=0$ and
$$\dot{Z} = AX_1 + \dot{U} - AX_2 - \dot{U} = A(X_1 - X_2) = AZ$$
But then we have that $Z = Z(t)$ is a solution to the homogeneous equation with initial condition $Z(0)=0$ so $Z = 0$ and then $X_1 = X_2$ which proves uniqueness. Now, with existence, we have just to find a solution than it will prove that solutions exists and are unique.
Solution
Let define the Dyson Series, given an matrix $A = A(t)$, $I$ is $n\times n$ identity matrix,
$$D(t,s):= I + \sum_{n=1}^{\infty}\int_s^{t}\int_s^{t_1}\dots\int_s^{t_{n-1}}A(t_1)\dots A(t_n)\mathrm dt_n \dots \mathrm dt_1$$
Then I affirm (now without a proof) that the solution to the equation $\dot{X} = A(t)X(t) + F(t)$ with initial condition $X(0) = X_0$ is
$$X(t) = D(t)X_0 + \int_0^tD(t,s)F(s)\mathrm ds$$
If $A$ is a constant matrix as in your problem we have that (try to prove that ) $D(t,s) = e^{A(t-s)}$. So the solution to your problem will be, given an initial condition $X_0 = X(0)$
$$X(t) = e^{At}X_0+\int_0^te^{A(t-s)}U(s)\mathrm ds$$
So, for the complete solution you have to find $e^A$. This is not a very symple task. Just because of the element $a_{21}$. If you make an $n$ product of the matrix $A$ you'll have
$$A^n = \left(\begin{matrix}a_{11}^n & 0 \\ f^n(a_{11},a_{21},a_{22}) & a_{22}^n\end{matrix}\right)$$
Where
$$f^1(a_{11}, a_{21}, a_{22}) = a_{21}$$
$$f^2(a_{11}, a_{21}, a_{22}) = a_{21}(a_{11}+a_{22})$$
$$f^3(a_{11}, a_{21}, a_{22}) = a_{21}(a_{11}^2+a_{22}(a_{11}+a_{22}))$$
