I will be going to set $\epsilon_0=\mu_0=1$. Now the Maxwell equations are:
$$\nabla . \textbf{E} = \rho$$
$$\nabla . \textbf{B} = 0$$
$$\nabla \times \textbf{E} = -\frac{\partial\textbf{B}}{\partial t}$$
$$\nabla \times \textbf{B} = \left(\textbf{J} + \frac{\partial\textbf{E}}{\partial t}\right)$$
and we have the identity
\begin{equation}
\nabla \times (\nabla \times\textbf{V}) = \nabla(\nabla.\textbf{V}) - \nabla^2\textbf{V}
\end{equation}
Now proceeding in the same fashion as in the vacuum
$$\nabla \times (\nabla \times \textbf{E}) = \nabla \times -\frac{\partial\textbf{B}}{\partial t} = -\frac{\partial}{\partial t}(\nabla \times \textbf{B})=-\frac{\partial}{\partial t}\left(\textbf{J} + \frac{\partial\textbf{E}}{\partial t}\right)$$
while the LHS becomes:
$$\nabla(\nabla.\textbf{E}) - \nabla^2\textbf{E}=\nabla(\rho) - \nabla^2\textbf{E}$$
Rearranging RHS and LHS we get
$$\nabla^2\textbf{E}-\frac{\partial^2\textbf{E}}{\partial t^2}=\nabla\rho +\frac{\partial}{\partial t}\textbf{J}$$
In simpler terms $$\Box\textbf{E}=\textbf{C}$$
where $$\textbf{C}=\nabla\rho +\frac{\partial}{\partial t}\textbf{J}$$
Now moving to the case of $\textbf{B}$
$$\nabla \times (\nabla \times \textbf{B})=\nabla \times\left(\textbf{J} + \frac{\partial\textbf{E}}{\partial t}\right)= \nabla \times\textbf{J} + \frac{\partial}{\partial t}(\nabla\times\textbf{E})=\nabla \times\textbf{J} -\frac{\partial^2\textbf{E}}{\partial t^2}$$
as for LHS we have
$$\nabla(\nabla.\textbf{B}) - \nabla^2\textbf{B}=\nabla(0) - \nabla^2\textbf{B}$$
Rearranging RHS and LHS we get
$$\nabla^2\textbf{B}-\frac{\partial^2\textbf{B}}{\partial t^2}=-\nabla \times\textbf{J} $$
in simpler terms $$\Box \textbf{B}=\textbf{F}$$
where $$\textbf{F}=-\nabla \times\textbf{J}$$
So putting sources has ultimately lead to what we call inhomogeneous wave equation which is simply $$\Box f(t,\vec{x})=h(t,\vec{x})$$ same thing as in case of Laplacian and Poisson equation in chapter 3.
Bonus material(I will make an assumption of tensors):
Maxwell equations are Lorentz covariant equation (that's how they contributed in Einstein triumph of special relativity), even when they were discovered in the era of Newtonian mechanics. Lorentz covariance another term to say a given physical quantity obey transformation law of different inertial reference frames in special relativity.
You might have also noticed how messy it becomes taking curl and div each time in above calculation and you will see it when you compare equation of chapters 10 and 12 of Griffiths book relating $\vec{J},\rho, A_\mu$. I will provide a rough sketch of the above calculation in the light of SR.
We define a quantity called 4-vector a generalization of vectors in 4 dimensions of Minkowski space $$A_{\mu}=(V, A_x, A_y, A_z)$$
$$J_{\mu}=(\rho, J_x, J_y, J_z)$$
Define a quantity called electromagnetic strength tensor $$F_{\mu\nu}=\partial_\mu A_\nu-\partial_\nu A_\mu$$
Maxwell equation can be recast as $$\partial^\nu F_{\mu\nu}=J_\mu$$ and $$\partial_{[\mu} F_{\nu\lambda]}=0$$
Leave the second equation aside (it's actually a tautology) let's focus on first equation expanding it in terms of $A_\mu$
$$\partial^\nu(\partial_\mu A_\nu-\partial_\nu A_\mu)=J_\mu$$
$$\partial^\nu(\partial_\mu A_\nu)-\partial^\nu(\partial_\nu A_\mu)=J_\mu$$ rearanging the terms we have
$$\partial_\mu(\partial^\nu A_\nu)-(\partial^\nu\partial_\nu) A_\mu=J_\mu$$
Now we use the Lorentz gauge and set $\partial^\nu A_\nu$ so ultimately we're left with
$$-(\partial^\nu\partial_\nu) A_\mu=J_\mu$$ which is nothing but $$-\Box A_\mu = J_\mu$$
which is simply the wave equation of different potentials under the presence of different sources you can recover $\vec{E}$, $\vec{B}$ from $A_\mu$. You might not have followed anything from this bonus material if it's your first encounter of,4-vectors, tensors, Einstein summation, Gauge transformation/freedom. What I actually wanted to show you was bear for the time being in the complicated mess of calculation and when you're done with chapter 12 of Griffith you'll have a different outlook at Electrodynamics as a whole.
Best Answer
There's nothing wrong with the first order wave equation mathematically, but it's just a little boring. If you want to use this equation to describe waves, it basically amounts to having a 1d solid with speed of sound $v$ for left moving waves (say) and speed of sound $0$ for right moving waves. It wouldn't surprise me if such a thing could be constructed (you would have to introduce some external fields to break time reversal invariance) but it is a very special system that we are not generically interested in.
Let's take the Fourier transforms of both equations to get the dispersion relationships. The normal second order equation gives \begin{equation} \omega^2=v^2 k^2 \end{equation} So for each frequency $\omega$ there are two allowed values of $k$, corresponding to right and left moving waves. Note that if we generalize the second order equation to include more spatial directions, there would be an infinite number of allowed $k$ values.
The first order equation meanwhile always has one allowed solution for a given frequency \begin{equation} \omega=v k \end{equation} So we get either right moving or left moving waves but not both. This restricts the allowed behavior, you can't have standing waves for example. If I try to generalize to higher dinensions, this equation picks out a single allowed $k$ for each frequency, so waves will only propagate along one very special direction.
Physically this is not what we would normally call a wave bc I only need one initial condition, not two. Usually dynamical systems can only be evolved given their initial position and velocity, but the first order equation needs only the initial position. (Or if you like, your equation is not a Hamiltonian system bc the phase space is odd dimensional).
Last but not least the first order equation necessarily picks out a preferred frame. By doing a boost I can change the sign of $v$, Thus the equation is not a good starting point for dealing with relativistic waves, which is one major application for the wave equation. (Of course you can have waves in materials that do pick out a preferred frame, and that is fine, but there you run into the problems above that you are looking at something with a preferred direction of motion as well).
(The Dirac equation gets around this by using spinor reps of the Lorentz group, but from your question I am supposing $f$ is a scalar).
Edit: rereading your question I see you want to have the $\pm$. Then you aren't looking at solutions of one single equation, you are looking at solutions to two equations and saying both are allowed. This is a little ugly for a few reasons. First,philosophically there should be one single equation for any system. Second, super positions don't solve either first order equation separately but do solve the second order equation. Third, the analogue of your idea for more spatial dimensions is to have an infinite set of first order equations, one for each direction.
On the other hand there is a way to rewrite the second order equation as two first order equations in a way which generalizes to any dimension, this is the way of the Hamiltonian and it is indeed a very useful thing to do in many situations.