Let us rephrase the question as follows.
Is there a Lagrangian density ${\cal L}$ that (1) implements the global $SO(2)$ symmetry manifestly, and (2) whose Euler-Lagrange equations yield Maxwell's equations in vacuum where there are no sources?
The answer is Yes. Below we will show this. Let the speed of light be $c=1$ from now on.
1) Field variables. The model has $2\times 3=6$ gauge potential fields ${\cal A}^a_i(\vec{x},t)$. Here $i=1,2,3$ are three spatial directions, and $a=1,2$ is an internal $SO(2)$ index. The gauge potential transforms
$${\cal A}^a_i\to \sum_{b=1}^2M^a{}_b {\cal A}^b_i $$
in the 2-dimensional fundamental representation of $SO(2)$, where
$$M^a{}_b =\left[\begin{array}{cc} \cos(\theta) & \sin(\theta) \\ -\sin(\theta) &\cos(\theta)\end{array}\right] \in SO(2), \qquad \sum_{b,c=1}^2 (M^{t})_a{}^b g_{bc} M^c{}_d = g_{ad}, \qquad g_{ab}\equiv \delta_{ab}.$$
The magnetic field $\vec{B}$ and electric field $\vec{E}$ are given by the curl of the gauge potential $\vec{\cal A}^a$,
$$\vec{\cal B}^a := \vec{\nabla} \times \vec{\cal A}^a, \qquad a=1,2, $$
where
$$ \vec{B}\equiv\vec{\cal B}^1 \qquad \mathrm{and} \qquad \vec{E}\equiv \vec{\cal B}^2.$$
It is easy to check that
$${\cal B}_i^a\to \sum_{b=1}^2 M^a{}_b {\cal B}_i^b$$
implements the sought-for $SO(2)$ transformation on the magnetic and electric fields $\vec{B}$ and $\vec{E}$. The two scalar-valued Maxwell equations
$$\vec{\nabla}\cdot\vec{B}=0 \qquad \mathrm{and} \qquad \vec{\nabla}\cdot\vec{E}=0 $$
are identically satisfied, because
$$\vec{\nabla}\cdot\vec{\cal B}^a=\vec{\nabla}\cdot(\vec{\nabla}\times\vec{\cal A}^a)=0, \qquad a=1,2.$$
2) Lagrangian density. The Lagrangian density ${\cal L}$ is
$$
{\cal L} = \frac{1}{2}\sum_{a,b=1}^2 \vec{\cal B}^a \cdot\left( \epsilon_{ab} \frac{\partial \vec{\cal A}^b }{\partial t} - g_{ab} \vec{\cal B}^b\right),
$$
where $\epsilon_{ab} = - \epsilon_{ba}$ is the Levi-Civita tensor in 2 dimensions (with $\epsilon_{12} = 1$). It is easy to check that the Lagrangian density ${\cal L}$ is manifestly invariant under global $SO(2)$ transformations, because both $\epsilon_{ab}$ and $g_{ab}\equiv \delta_{ab}$ are invariant tensors for $SO(2)$. The action is by definition
$$ S[{\cal A}]=\int d^4x \ {\cal L}. $$
Extremizing the action with respect to $2\times 3=6$ gauge potential fields ${\cal A}^a_i$ produces $6$ Euler-Lagrange equations. In detail, an arbitrary infinitesimal variation ${\cal A}^a_i\to {\cal A}^a_i+\delta{\cal A}^a_i$ induces a change $\delta{\cal L}$ in the Lagrangian density
$$
\delta{\cal L} ~\sim~ \sum_{a,b=1}^2\delta \vec{\cal B}^a \cdot\left( \epsilon_{ab} \frac{\partial \vec{\cal A}^b }{\partial t} - g_{ab} \vec{\cal B}^b\right)
~\sim~ \sum_{a,b=1}^2\delta \vec{\cal A}^a\cdot \left( \epsilon_{ab} \frac{\partial \vec{\cal B}^b }{\partial t} - g_{ab} \vec{\nabla} \times \vec{\cal B}^b\right),
$$
where the "$\sim$" sign means equality modulo divergence terms. The variation $\delta S=0$ of the action vanishes iff the last parenthesis is zero. This yield precisely the two vector-valued Maxwell equations
$$ \frac{\partial \vec{B}}{\partial t} + \vec{\nabla} \times \vec{E}=\vec{0} \qquad \mathrm{and} \qquad \frac{\partial \vec{E}}{\partial t} = \vec{\nabla} \times \vec{B}.
$$
References:
1) S. Deser and C. Teitelboim, "Duality Transformations Of Abelian And Nonabelian Gauge Fields", Phys.Rev.D 13 (1976) 1592.
2) C. Bunster and M. Henneaux, "Can (Electric-Magnetic) Duality Be Gauged?", Phys.Rev.D83 (2011) 045031, arXiv:1011.5889.
3) S. Deser, "No local Maxwell duality invariance", Class.Quant.Grav.28 (2011) 085009, arXiv:1012.5109.
Best Answer
In a (classical) lagrangian field theory, the configuration space $\mathcal C$ of the system is a space of field configurations. A field configuration (or just "field" for short) is usually taken to be a function $\phi:\mathcal M\to T$ where $M$ is a manifold and $T$ is some set, often a manifold or vector space or both, called the target space of the field. The configuration space $\mathcal C$ is then taken to be some sufficiently smooth (when a notion of smoothness can be defined) subset of the set of all possible fields. The lagrangian is then a function $L:\mathcal C\to\mathcal C$; \begin{align} \phi\mapsto L[\phi] \end{align} Namely, the lagrangian maps a particular field configuration, to another field configuration. Often, one considers a field theory for which the lagrangian can be written as local local density, but this is not strictly speaking necessary. The action of the theory can then be defined as the integral of $L[\phi]$ over $M$; \begin{align} S[\phi] = \int_M \, d^Dx\,L[\phi](x). \end{align} Note. My terminology and notation here are a bit non-standard in some contexts. For example, in relativistic physics (field theory) the Lagrangian will usually map a field configuration $\phi$ to a function $L[\phi]$ of time, and then this function of time will be integrated to yield the action. It's not hard in practice to translate between conventions.
One can then define what it means for the action to possess symmetry. In particular, given a mapping $F:\mathcal C\to \mathcal C$ of the manifold on which field configurations are defined to itself, one says that the action is invariant under $F$ provided \begin{align} S[F(\phi)] = S[\phi] \end{align} for all $\phi\in\mathcal C$. For "continuous" transformations, one can also define notions of symmetry that don't involve full invariance, but let's keep the discussion simple at this point.
Example. A common toy theory considered as the first example in most relativistic field theory texts is that of a single, free, real Lorentz scalar defined on Minkowski space (I'll use metric signature $+ - - -$). In this case, we have \begin{align} M = \mathbb R^{3,1}, \qquad T = \mathbb R \end{align} and $\mathcal C$ is a space of sufficiently smooth functions $\phi:\mathbb R^{3,1}\to\mathbb R$ that satisfy certain desired boundary conditions. The Lagrangian of such a theory is \begin{align} L[\phi](x) = \mathscr L(\phi(x), \partial_0\phi(x), \dots \partial_3\phi(x)) \end{align} where $\mathscr L$ is the Lagrangian density defined by \begin{align} \mathscr L(\phi, \partial_0\phi, \dots, \partial_3\phi) &= \frac{1}{2}\partial_\mu\phi\partial^\mu\phi - \frac{1}{2}m^2\phi^2. \end{align} Given any Lorentz transformation $\Lambda$, one can define a transformation $F_\Lambda:\mathcal C\to\mathcal C$ as follows: \begin{align} F_\Lambda(\phi)(x) = \phi(\Lambda^{-1} x). \end{align} A short computation then shows that the Lagrangian is a Lorentz scalar under this transformation, namely \begin{align} L[F_\Lambda(\phi)](x) = L[\phi](\Lambda^{-1}x). \end{align} In fact, this is essentially done for you on page 36 of Peskin. It follows from this that the action is invariant under $F$; \begin{align} S[F_\Lambda(\phi)] = \int_{\mathbb R^{3,1}} d^4x\, L[\phi](\Lambda^{-1}x) = \int_{\mathbb R^{3,1}} d^4x\, L[\phi](x) = S[\phi]. \end{align} since the measure $d^4 x$ is Lorentz-invariant. Notice, in particular, that the fact that the Lagrangian transformed as a Lorentz scalar (namely precise in the same way as the scalar field $\phi$ was defined to transform), immediately led to invariance of the action.
Furthermore, suppose that $\phi$ is a field configuration that leads to stationary action, then we can also show that a Lorentz-transformed field leads to a stationary action using Lorentz invariance of the action. To see this, recall that the variational derivative in the direction of a field configuration $\eta$ is defined as follows: \begin{align} \delta_\eta S[\phi] = \frac{d}{d\epsilon}S[\phi+\epsilon\eta]\Big|_{\epsilon=0} \end{align} Now, suppose that $\phi$ is a stationary point of the action, namely that $\delta_\eta S[\phi] = 0$ for all admissible $\eta$, then for all such $\eta$ we have \begin{align} \delta_{F_\Lambda(\eta)} S[F_\Lambda(\phi)] &= \frac{d}{d\epsilon} S[F_\Lambda(\phi) + \epsilon F_\Lambda(\eta)]\Big|_{\epsilon = 0} \\ &= \frac{d}{d\epsilon} S[F_\Lambda(\phi+\epsilon_\eta)]\Big|_{\epsilon = 0} \\ &= \frac{d}{d\epsilon} S[\phi+\epsilon_\eta]\Big|_{\epsilon = 0} \\ &= 0 \end{align} Now set $\eta = F_\Lambda^{-1}(\xi)$, then the computation we just performed shows that \begin{align} \delta_{\xi} S[F_\Lambda(\phi)] =0. \end{align} for all admissible field configurations $\xi$. In other words, the Lorentz transformed scalar is also a stationary point of the action. Notice that this demonstration holds for any Lorentz transformation, not just boosts.
Addendum.
As pointed out in the comments, the argument at the end about variational derivatives hinges on linearity of $F_\Lambda$. This can be demonstrated as follows: \begin{align} F_\Lambda(a\phi+b\psi)(x) &= (a\phi+b\psi)(\Lambda^{-1}x) \\ &= a\phi(\Lambda^{-1}x) + b\psi(\Lambda^{-1}x) \\ &= aF_\Lambda(\phi)(x) + bF_\Lambda(\psi)(x). \end{align}
Let me make some remarks about the mapping $F:\mathcal C\to \mathcal C$; a symmetry of the action. If there exists a mapping $f_T:T\to T$ on the target space that induces this mapping, namely \begin{align} F(\phi)(x) = f_T(\phi(x)), \end{align} then $F$ is called an internal symmetry. On the other hand, if there is a mapping $f_M:M\to M$ on the base manifold $M$ that induces this mapping, namely \begin{align} F(\phi)(x) = \phi(f_M(x)), \end{align} then $F$ is called a base manifold symmetry (or more commonly a spacetime symmetry since in the context of relativistic field theory, the base manifold is a spacetime like Minkowski space.)
Furthermore, the mapping $F:\mathcal C \to\mathcal C$ on the field configuration space is often, as in the scalar field example, a group action of some group $G$ on $\mathcal C$. This means that to each $g\in G$, we associate a mapping $F_g:\mathcal C\to \mathcal C$ such that the mapping $g\mapsto F_g$ is a homomorphism of the group $G$. In practice, the group $G$ is sometimes a group of symmetries that naturally acts on the base manifold, and sometimes $G$ a group of symmetries that naturally acts on the target space (or even both when the $M=T$). In any event, this group action is usually obtained by composing a target space group action $(f_T)_g:T\to T$ with a base manifold group action $(f_M)_g:M\to M$. More explicitly, for each $g\in G$, we can define mappings $(F_T)_g:\mathcal C\to \mathcal C$ and $(F_M)_g:\mathcal C\to\mathcal C$ as follows: \begin{align} (F_T)_g(\phi)(x) = (f_T)_g(\phi(x)), \qquad (F_M)_g(\phi)(x) = \phi((f_M)_g(x)) \end{align} and then the full group action $F_g:\mathcal C\to\mathcal C$ is defined by the composition of these two; \begin{align} F_g = (F_T)_g\circ (F_M)_g \end{align} or more explicitly \begin{align} F_g(\phi)(x) = (f_T)_g(\phi((f_M)_g(x))). \end{align} Now, this is all a bit abstract, so let's write out what all of these objects would be for the scalar field example: \begin{align} G &= \mathrm{SO}(3,1) \\ g &= \Lambda\\ (f_T)_g(\phi(x)) &= \phi(x) \\ (f_M)_g(x) &= \Lambda^{-1}x \\ (F_T)_g(\phi)(x) &= \phi(x) \\ (F_M)_g(\phi)(x) &= \phi(\Lambda^{-1}x) \\ F_g(\phi)(x) &= \phi(\Lambda^{-1}x) \end{align} Notice that $f_T$ is simply the identity mapping on the target space. This is precisely what we mean when we say that the scalar field is a Lorentz scalar. On the other hand, for a Lorentz vector, the target space itself would be Minkowski space $\mathbb R^{3,1}$, and the target space group action would be \begin{align} (f_T)_\Lambda(A(x)) = \Lambda A(x), \end{align} namely, there is an internal symmetry in which the vector indices on the field transform non-trivially. In components, which is how you'll see this written in Peskin for example, the right hand side would be written as $\Lambda^\mu_{\phantom\mu\nu} A^\mu(x)$.