I) Let us just for fun generalize OP's question to $n$ spacetime dimensions, and check how the counting of eqs. and degrees of freedom (d.o.f.) work out in this general setting. We shall use Lubos Motl's answer as a template for this part. Also we shall use a special relativistic $(-,+,\ldots,+)$ notation with $c=1$, where $\mu,\nu\in\{0,\ldots,n-1\}$ denote spacetime indices, while $i,j \in\{1,\ldots,n-1\}$ denote spatial indices. Maxwell eqs. are the following.
Source-free Bianchi identities:
$${\rm d}F~=~0 \qquad\qquad \Leftrightarrow \qquad\qquad
\sum_{\rm cycl.~\mu,\nu,\lambda} d_{\lambda} F_{\mu\nu} ~=~0, \qquad\qquad
F~:=~\frac{1}{2} F_{\mu\nu}~ {\rm d}x^{\mu} \wedge {\rm d}x^{\nu}.$$
Here
$$\left(\begin{array}{c} n \cr 3\end{array}\right) {\rm~Bianchi~identities} ~=~ \left(\begin{array}{c} n-1 \cr 3\end{array}\right) {\rm~constraints}~+~ \left(\begin{array}{c} n-1 \cr 2\end{array}\right) {\rm~dynamical~eqs.} $$
$$~=~ ({\rm No~magnetic~monopole~eqs.})~+~ ({\rm Faraday's~law}). $$
Maxwell eqs. with source terms:
$$ d_{\mu}F^{\mu\nu}~=~-j^{\nu} .$$
Here
$$n {\rm~source~eqs.}~=~1 {\rm~constraint} ~+~ (n-1) {\rm~dynamical~eqs.}$$
$$~=~({\rm Gauss'~law}) ~+~ ({\rm Ampere's~law~with~displacement~term}).$$
We have used the terminology that a dynamical eq. contains time derivatives, while a constraint does not. So the number of dynamical eqs. is
$$ \left(\begin{array}{c} n-1 \cr 2\end{array}\right)~+~(n-1)~=~ \left(\begin{array}{c} n \cr 2\end{array}\right),$$
which precisely matches
$${\rm the~number~} \left(\begin{array}{c} n \cr 2\end{array}\right) {\rm~of~}
F_{\mu\nu} {\rm~fields}$$
$$~=~\left(\begin{array}{c} n-1 \cr 2\end{array}\right){~\rm magnetic~fields~} F_{ij} ~+~(n-1) {\rm~electric~fields~}F_{i0} .$$
Maxwell eqs. with source terms imply the continuity eq.
$$ d_{\nu}j^{\nu} ~=~-d_{\nu}d_{\mu}F^{\mu\nu}~=~0,\qquad\qquad F^{\mu\nu}~=~-F^{\nu\mu},$$
so one must demand that the background sources $j^{\nu}$ obey the continuity eq.
For consistency, the time derivative of each of the constraints should vanish. In the case of the no-magnetic-monopole-eqs., this follows from Faraday's law. In the case of Gauss' law, this follows from the modified Ampere's law and the continuity eq.
II) The previous section (I) made the counting in terms of the $\left(\begin{array}{c} n \cr 2\end{array}\right)$ field strengths $F_{\mu\nu}$. In terms of the $n$ gauge potentials $A_{\mu}$, the counting goes as follows. The Bianchi identities are now trivially satisfied,
$$F~=~{\rm d}A\qquad\qquad A~:=~A_{\mu}~ {\rm d}x^{\mu}. $$
There are still the $n$ Maxwell eqs. with source terms
$$ (\Box\delta^{\mu}_{\nu}-d^{\mu}d_{\nu})A^{\nu}~=~-j^{\mu} , \qquad\qquad \Box~:=~d_{\mu}d^{\mu}. $$
There is a single gauge d.o.f. because of gauge symmetry $A \to A + {\rm d}\Lambda$ and $F \to F$. If one gauge-fixes using the Lorenz gauge condition
$$d_{\mu}A^{\mu}~=~0, $$
the Maxwell eqs. become $n$ decoupled wave equations
$$ \Box A^{\mu}(x)~=~-j^{\mu}(x). $$
By a spatial Fourier transformation, these become decoupled linear second-order ODEs with constant coefficients,
$$ (d^2_t+\vec{k}^2) \hat{A}^{\mu}(t;\vec{k})~=~\hat{j}^{\mu}(t;\vec{k}) , $$
which, starting from some initial time $t_0$, may be solved for all times $t$, cf. OP's question. [One should check that the solution
$$\hat{A}^{\mu}(t;\vec{k})
~=~\int {\rm d} t^{\prime} ~G(t-t^{\prime};\vec{k})~\hat{j}^{\mu}(t^{\prime};\vec{k}), \qquad\qquad
(d^2_t+\vec{k}^2)G(t-t^{\prime};\vec{k})~=~\delta(t-t^{\prime}),$$
satisfies the Lorenz gauge condition. This follows from the continuity eq.]
III) It is interesting to derive the complete solution $\tilde{A}^{\mu}(k)$ in $k^{\nu}$-momentum space without gauge-fixing. The Fourier-transformed Maxwell eqs. read
$$M^{\mu}{}_{\nu}~\tilde{A}^{\nu}(k)~=~\tilde{j}^{\mu}(k), \qquad\qquad M^{\mu}{}_{\nu}~:=~k^2\delta^{\mu}_{\nu} -k^{\mu}k_{\nu}. $$
To proceed one must analyze the matrix $M^{\mu}{}_{\nu}$ for fixed $k^{\lambda}$. There are three cases.
Constant mode $k^{\mu}=0$. Then the matrix $M^{\mu}{}_{\nu}=0$ vanishes identically. Maxwell eqs. are only possible to satisfy if $\tilde{j}^{\mu}(k=0)=0$ is zero. The gauge potential $\tilde{A}_{\mu}(k=0)$ is not restricted at all by Maxwell eqs., i.e., there is a full $n$-parameter solution.
Massive case $k^2\neq 0$. The matrix $M^{\mu}{}_{\nu}$ is diagonalizable with eigenvalue $k^2$ (with multiplicity $n-1$), and eigenvalue $0$ (with multiplicity $1$). The latter corresponds to a pure gauge mode $\tilde{A}^{\mu}~\propto~k^{\mu}$. The complete solution is a $1$-parameter solution of the form
$$\tilde{A}^{\mu}(k)
~=~\frac{\tilde{j}^{\mu}(k)}{k^2}~+~ik^{\mu}\tilde{\Lambda}(k).$$
Apart from the source term, this is pure gauge.
Massless case $k^2=0$ and $k^{\mu}\neq 0$. The matrix $M^{\mu}{}_{\nu}$ is not diagonalizable. There is only eigenvalue $0$ (with multiplicity $n-1$).
Maxwell eqs. are only possible to satisfy if the source $\tilde{j}^{\mu}(k)=\tilde{f}(k)k^{\mu}$ is proportional to $k^{\mu}$ with some proportionality factor $\tilde{f}(k)$. In that case Maxwell eqs. become
$$ -k_{\mu}\tilde{A}^{\mu}(k)~=~\tilde{f}(k). $$
Let us introduce an $\eta$-dual vector$^1$
$$k^{\mu}_{\eta}~:=~(-k^0,\vec{k})\qquad {\rm for}\qquad
k^{\mu}~=~(k^0,\vec{k}).$$
Note that
$$k_{\mu}~k^{\mu}_{\eta}~=~(k^0)^2+\vec{k}^2$$
is just the Euclidean distance square in $k^{\mu}$-momentum space. The complete solution is an $(n-1)$-parameter solution of the form
$$\tilde{A}^{\mu}(k)
~=~-\frac{k^{\mu}_{\eta}}{k_{\nu}~k^{\nu}_{\eta}}\tilde{f}(k)
~+~ik^{\mu}\tilde{\Lambda}(k)~+~\tilde{A}^{\mu}_{T}(k).$$
The term proportional to $k_{\mu}$ is pure gauge. Here $\tilde{A}^{\mu}_{T}(k)$ denote $n-2$ transversal modes,
$$k_{\mu}~\tilde{A}^{\mu}_{T}(k)~=~0, \qquad\qquad
k_{\mu}^{\eta}~\tilde{A}^{\mu}_{T}(k)~=~0. $$
The $n-2$ transversal modes $\tilde{A}^{\mu}_{T}$ are the only propagating physical d.o.f. (electromagnetic waves, photon field).
--
$^1$ Longitudinal and timelike polarizations are in the massless case proportional to $k^{\mu}\pm k^{\mu}_{\eta}$, respectively.
