I think it's important to emphasize a point which was never emphasized to me when I was taking my courses. When a theory such as electromagnetism can be formulated in (roughly) equivalent ways at varying levels of sophistication, the more sophisticated forms should make your life easier, not harder.
Now, generally speaking this benefit only comes when you actually understand the added sophistication from both a physical and mathematical point of view, so for a time it may seem more complicated. But when that understanding has been achieved, there will be a moment where you suddenly see through the symbols and recognize everything as simpler than you first thought. The higher mathematics is there to serve you, not the other way around.
To put this in practice with electromagnetism, start from a relativistic standpoint and imagine that there was a force which is a linear function of the 4-velocity - that is, something of the form
$$\dot{\mathbf p} = qF \mathbf u$$
where $\mathbf p = m\mathbf u$ is the 4-momentum, $\mathbf u$ is the 4-velocity, and $q$ is some coupling constant. What are the properties of this $F$, which must be expressed as a matrix of some kind? Well, the first thing to point out is that because $\mathbf p\cdot \mathbf p = m^2c^4$ is constant, it follows that $\mathbf p \cdot \dot{\mathbf p} = 0$. But this requires that $$\mathbf u \cdot F\mathbf u = u^\mu\big(g_{\mu \alpha} F^\alpha_{\ \ \nu} \big)u^\nu \equiv F_{\mu\nu} u^\mu u^\nu =0$$ for any possible (timelike) 4-vector $\mathbf u$. The simplest way to satisfy this is to require that $F_{\mu\nu}$ be antisymmetric.
In any specific frame of reference, such a matrix acting on the 4-velocity can be broken up into a piece which depends on the ordinary 3-velocity $\vec v$ (in that frame) and a piece which does not. It's a straightforward exercise to show that
$$qF\mathbf u = q\gamma \pmatrix{\sum_{i=1}^3 F^0_{\ \ i} v^i \\ F^1_{\ \ 0} + \sum_{i=1}^3 F^1_{\ \ i} v^i \\ F^2_{\ \ 0} + \sum_{i=1}^3 F^2_{\ \ i} v^i \\ F^3_{\ \ 0} + \sum_{i=1}^3 F^3_{\ \ i} v^i }$$
This suggests that, for organizational purposes, we might define
$$F^i_{\ \ 0}\equiv E^i \qquad F^i_{\ \ j} = \mathcal B^i_{\ \ j}$$
where $\vec E$ is a 3-vector and $\mathcal B$ is a 3$\times$3 matrix with the property that $\mathcal B^i_{\ \ j} = -\mathcal B^j_{\ \ i}$ for each $i,j$. With this notation being defined, we have
$$F\mathbf u = q\gamma \pmatrix{\vec E \cdot \vec v \\ \vec E + \mathcal B\vec v}$$
Finally, observe that in 3 dimensions, the action of an antisymmetric matrix on a vector can be re-expressed as as a cross product: $$\pmatrix{0 & \mathcal B^1_{\ \ 2} & \mathcal B^1_{\ \ 3} \\ -\mathcal B^1_{\ \ 2} & 0 & \mathcal B^2_{\ \ 3}\\ -\mathcal B^1_{\ \ 3} & -\mathcal B^2_{\ \ 3} & 0}\pmatrix{v^1\\v^2\\v^3} = \pmatrix{\mathcal B^1_{\ \ 2}v^2 + \mathcal B^1_{\ \ 3} v^3 \\ - \mathcal B^2_{\ \ 1} v^1 + \mathcal B^2_{\ \ 3}v^3\\ -\mathcal B^1_{\ \ 3} v^1 - \mathcal B^2_{\ \ 3} v^2 } = \vec v \times \vec B $$ where we have defined the 3-vector $\vec B \equiv (-\mathcal B^2_{\ \ 3}, \mathcal B^1_{\ \ 3}, -\mathcal B^1_{\ \ 2})$. So finally we have obtained
$$F\mathbf u = \gamma \pmatrix{q\vec E \cdot \vec v\\ q\big(\vec E + \vec v\times \vec B\big)}$$
The main lessons from this go as follows.
- The fact that the force is linear in $\mathbf u$ and preserves the rest mass of the particle is sufficient to show that in any specific reference frame, it can be decomposed into timelike and spacelike parts - or put differently, into parts proportional to $\vec v$ and parts which are not.
- This splitting is obviously frame-dependent, because $\vec v$ changes from frame to frame.
- The $\vec v$-dependent part is the action of an antisymmetric matrix $\mathcal B$ on $\vec v$. The fact that this can be expressed as the cross product of $\vec v$ with a vector $\vec B$ is an accident of 3-dimensions, essentially deriving from the fact that an $n\times n$ antisymmetric matrix has $n(n-1)/2$ independent components.
- If we want Lorentz-invariant quantities, we should obtain them at the level of $F$. Examples include $F^2 \equiv F_{\mu\nu}F^{\mu \nu} = 2(B^2-E^2/c^2)$ and $\epsilon^{\mu\nu\alpha\beta}F_{\mu\nu}F_{\alpha\beta} \propto \vec E \cdot \vec B$.
If we wanted to go further and endow this $F$ with dynamics of its own, we could start to think about the simplest, Lorentz-covariant differential equations we could write down. In what may or may not be surprising, the Maxwell equations take the rather straightforward form $$\nabla_\mu F^{\mu\nu} = \mu_0 J^\nu \qquad \epsilon^{\mu \alpha\beta\gamma}\nabla_\alpha F_{\beta\gamma} = 0$$ where indices have been freely lowered and raised with the metric. The latter can be re-interpreted by defining the dual tensor $$G^{\alpha\mu} = \frac{1}{2}\epsilon^{\alpha \mu \beta \gamma}F_{\beta \gamma}$$ in which case the second equation becomes $\nabla_\mu G^{\mu\nu} = 0$. These equations can be discussed and motivated at various levels as well, but that's beyond the scope of this answer.
