Is there a combined Lagrangian that gives both the Lorentz force and Maxwell equations via the Euler-Lagrange equations?
Not sure about the continuous charge/current density case, but for point particles, this can be done.
There is the standard textbook Lagrangian for $N$ charged point particles:
$$
L(A^\nu,\partial_\mu A^\nu,\{\mathbf r_a,\mathbf v_a\}_{a=1}^{N}) = \int_V -\frac{1}{4\mu_0}F^{\mu\nu}F_{\mu\nu} d^3 \mathbf x~~~+$$
$$- ~~\sum_a q_a \varphi(\mathbf r_a,t) + \sum_a q_a\mathbf v_a \cdot \mathbf A(\mathbf r_a,t) ~~~+
$$
$$
-~~~\sum_a \sqrt{1-v_a^2/c^2}m_a c^2,
$$
where $\varphi,\mathbf A, F$ refer to total EM field. Thus this is a function of $N$ positions, $N$ velocities, and a functional of the fields $\varphi,\mathbf A$.
This Lagrangian is in textbooks used to "derive" (e.g. in Landau&Lifshitz) both the Maxwell equations for total fields in presence of the current density $\sum_a q_a\mathbf v_a \delta(\mathbf x - \mathbf r_a)$ and charge density $\sum_a q_a\delta(\mathbf x - \mathbf r_a)$, and also to "derive" the equations of motion for all the particles, with each particle experiencing the Lorentz force $q_a\mathbf E(\mathbf r_a,t) + q_a\mathbf v_a\times \mathbf B(\mathbf r_a,t)$.
The problem with this is that when $\varphi,\mathbf A$ refer to total fields (due to all particles), expressions such as $\varphi(\mathbf r_a,t)$ and $\mathbf v_a \cdot \mathbf A(\mathbf r_a,t)$ are undefined, because the fields $\varphi,\mathbf A$ are singular at positions of the particles. The equations of motion appear to be correct, but are not, because they refer to undefined forces $q_a\mathbf E(\mathbf r_a,t) + q_a\mathbf v_a\times \mathbf B(\mathbf r_a,t)$, because fields $\mathbf E,\mathbf B$ too are undefined at the positions of the particles.
One could try to regard the field functions $\varphi,\mathbf A,\mathbf E,\mathbf B$ as "doctored", to have just the right values at positions of the particles to get the actual forces acting on the particles. But even if this was mathematically possible (the functions are necessarily discontinuous), the "doctored" values cannot be found from the Lagrangian function, the resulting equations of motion and initial conditions; thus the Lagrangian is not predictive of the resulting motion. Without knowing the actual motion of the particles, the completion values can be anything and produce arbitrary forces.
The Lagrangian function above is mathematically valid and predictive when $\varphi, \mathbf A, F$ refer to external fields not due to particles themselves, e.g. a free EM field, or a field of distant sources not considered part of the system. But then such Lagrangian can be physically valid only approximately, because it ignores interaction between the particles in the system. This is fine if the particles are far from each other and interact weakly, so the dominant force is the one due to the external field, but is not if there are strongly interacting particles.
If we want to have a valid Lagrangian (for a system of point particles) which takes into account mutual interactions of the particles, we have to formulate it in a way which makes all the quantities in it defined and which produces equations of motions with all values of forces there defined.
This is the case for the Frenkel theory, where 1) all particles produce their own EM field obeying the Maxwell equations with point sources; this field is singular at the particle producing it, and thus different from any other EM field produced by the other particles; 2) every particle experiences Lorentz force due to all EM fields except the one it itself produces; no self-interaction of particles is allowed.
Thus we have $N$ potentials $\varphi_b$, and $N$ vector potentials $\mathbf A_b$. The Lagrangian respecting 1), 2) is
$$
L(\{A_b^\nu,\partial_\mu A_b^\nu\}_{b=1}^{N},\{\mathbf r_a,\mathbf v_a\}_{a=1}^{N}) ~~~=
$$
$$
=~~~ \int_V -\sum_{a,b}'\frac{1}{2\mu_0}F_a^{\mu\nu}F_{b,\mu\nu} ~d^3 \mathbf x~~~+$$
$$- ~~\sum_{a,b}' q_a \varphi_b(\mathbf r_a,t) + \sum_{a,b}' q_a\mathbf v_a \cdot \mathbf A_b(\mathbf r_a,t) ~~~+
$$
$$
-~~~\sum_a \sqrt{1-v_a^2/c^2}~m_a c^2.
$$
The summation is over all pairs of indices $a,b$, where each pair is present only once; the prime symbol next to the summation symbol means the case $a=b$ is to be omitted.
The Euler-Lagrange equations for coordinates $\mathbf r_a$ produce the expected equations of motion, free of self-interaction:
$$
\frac{d}{dt}\bigg(\gamma_a m_a \mathbf v_a\bigg) = q_a \sum_{b}' \mathbf E_b(\mathbf r_a,t) + q_a\mathbf v_a \times \sum_{b}' \mathbf B_b(\mathbf r_a,t).
$$
where $\mathbf E_b$ is a short-hand for $-\nabla\varphi_b -\partial_t \mathbf A_b$, and $\mathbf B_b$ is a short-hand for $\nabla\times\mathbf A_b$, and they also produce all the Maxwell equations for every pair $\mathbf E_b,\mathbf B_b$.
The source-free Maxwell equations (the Gauss law for magnetic field and the Faraday law) follow from definition of the fields $\mathbf E_b,\mathbf B_b$ in terms of $\varphi_b,\mathbf A_b$, given above. This definition is motivated by the requirement that net force in the equations of motion of any particle should be the well-known Lorentz force due to electric and magnetic field due to all the other particles. Also, each four-current $j_a^\mu$ is locally conserved, because $q_a$ is assumed constant in time. The source Maxwell equations (the Gauss law for electric field and the Maxwell-Ampere law) for all pairs $\mathbf E_b,\mathbf B_b$ come out directly from the Euler-Lagrange equations for generalized coordinates $A_b^\mu$, when we express the derivatives of $A_b$ in terms of the above defined $\mathbf E_b$,$\mathbf B_b$.
See also the short paper
R. C. Stabler, A Possible Modification of Classical Electrodynamics, Physics Letters, 8, 3, (1964), p. 185-187. http://dx.doi.org/10.1016/S0031-9163(64)91989-4
and the original paper by Frenkel
J. Frenkel, Zur Elektrodynamik punktförmiger Elektronen, Zeits. f. Phys., 32, (1925), p. 518-534. http://dx.doi.org/10.1007/BF01331692