I would be interested in a detailed description of Maxwell's equations from a mathematical point of view, that is as first-order partial differential equations. Taking the equations in SI units in vacuum, $$ \begin{align} \nabla \cdot \mathbf{E}(\mathbf{r}, t) &= \frac{\rho(\mathbf{r}, t)}{\epsilon_0} \tag{1} \label{eq 1}\\ \nabla \cdot \mathbf{B}(\mathbf{r}, t) &= 0 \tag{2} \label{eq 2} \\ \frac{\partial \mathbf{B}(\mathbf{r}, t)}{\partial t} &= -\nabla \times \mathbf{E}(\mathbf{r}, t) \tag{3} \label{eq 3} \\ \frac{\partial \mathbf{E}(\mathbf{r}, t)}{\partial t} &= \frac{1}{\epsilon_0 \mu_0} \nabla \times \mathbf{B}(\mathbf{r}, t) - \frac{1}{\epsilon_0} \mathbf{j}(\mathbf{r}, t) \tag{4} \label{eq 4} \end{align} $$ where $\mathbf{r} = (x, y, z)$. As a first step, I will assume the functions $\rho(\mathbf{r}, t)$ and $\mathbf{j}(\mathbf{r}, t)$ to be prescribed, but arbitrary functions not affected by the electromagnetic field. Ignoring equations \eqref{eq 1} and \eqref{eq 2} for the moment, the remaining two equations form a system of six linear first-order coupled partial differential equations for six unknown functions $E_i(\mathbf{r}, t)$ and $B_i(\mathbf{r}, t)$ of four independent variables $x, y, z, t$.
First of all, is there a general theorem for partial differential equations that guarantees the existence of a solution for $\mathbf{E}$ and $\mathbf{B}$ for arbitrary $\mathbf{j}(\mathbf{r}, t)$? Secondly, if solutions for given $\mathbf{j}$ exist, what data is needed to uniquely determine a solution? In particular, is the specification of $\mathbf{E}(\mathbf{r}, t_0)$ and $\mathbf{B}(\mathbf{r}, t_0)$ for all $\mathbf{r}$ at some time $t_0$ enough?
Now add equations \eqref{eq 1} and \eqref{eq 2} . The system of equations now consists of eight linear first-order coupled partial differential equations for the six unknown functions. What about the existence and uniqueness of solutions now? Starting from \eqref{eq 3} and \eqref{eq 4}, it seems to me that adding \eqref{eq 1} and \eqref{eq 2} amounts to adding constraints on the possible solutions of \eqref{eq 3} and \eqref{eq 4}. Is there a guarantee that these constraints, for arbitrary $\rho(\mathbf{r}, t)$, are compatible with general solutions to \eqref{eq 3} and \eqref{eq 4}? Does adding these constraints reduce the number of functions which have to be solved for so that a reduced system of differential equations can be written down?
Essentially, I am interested in what a mathematician would have to say if I gave him these coupled equations for six unknown functions $u_{1,...,6}(\mathbf{r}, t)$ (the components of the fields), without any physical interpretation.
Ideally, I would be interested in answers to these questions, with as much mathematical rigour as possible, that only involve the fields $\mathbf{E}$ and $\mathbf{B}$. A description in terms of scalar or vector potentials would also be interesting, but only in addition to a treatment in terms of the fields themselves.
This is a partial answer. The answer to your second question ("Secondly ...") can be found on Page 647 of Courant and Hilbert, Methods of Mathematical Physics Volume II. The answer is yes, for a given current and initial values, solutions to (3) and (4) are unique. A proof is given there, that when the current is zero, solutions are unique. This answers your question because if you have two solutions with the same current and initial values, their difference satisfies the vacuum equations with zero initial values. That is what is proved to be zero in the book.