A point charge q appears at the origin at time 0, in an infinite vacuum. What's the electric field a billion miles away at time 0?
I suspect the real answer is 0. But it seems to me that Maxwell's first equation gradient(E) = (charge density) give us the answer q/(a billion miles) * (unit vector pointing away from the origin). Did I go wrong somewhere, or is this an example of Maxwell's equations not modelling nature perfectly?
The comments have already sketched the answer to the question, but I thought I might flesh it out with a bit more detail.
$$\nabla \cdot \vec E = \rho /\epsilon_0$$ $$\nabla \cdot \vec B = 0$$ $$\nabla \times \vec E = -\frac{\partial}{\partial t} \vec B$$ $$\nabla \times \vec B = \mu_0 \vec J + \epsilon_0 \mu_0 \frac{\partial }{\partial t} \vec E$$
First - why can't a single point charge appear out of nowhere? The answer is that the continuity equation is actually embedded in Maxwell's equation. Behold:
$$ \frac{\partial \rho}{\partial t} = \epsilon_0 \nabla \cdot \left(\frac{\partial }{\partial t} \vec E \right) = \epsilon_0 \nabla \cdot \left(\frac{\nabla \times \vec B}{\epsilon_0 \mu_0} - \frac{1}{\epsilon_0} \vec J\right) = -\nabla \cdot \vec J$$ so $$ \frac{\partial \rho}{\partial t} + \nabla \cdot \vec J = 0$$ where we've used the fact that the divergence of a curl is always zero. The meaning of this equation becomes more clear if we integrate it over an arbitrary volume $V$ (with boundary $\partial V$) and then apply the divergence theorem:
$$ \frac{d}{dt}Q_{in} = - \oint_{\partial V} \vec J \cdot d\vec S$$
In words, the rate of change of the total charge in any volume is equal to the negative of the current flux through the boundary of that volume. That is, any change in the total charge in any volume can be attributed to charge flowing through the boundary. Isolated charges cannot simply pop into existence.
Now to the actual content of your question. The best route to the answer involves introducing the scalar potential $\phi$ and the vector potential $\vec A$.
Because $\nabla \cdot \vec B = 0$ everywhere, we can write $\vec B = \nabla \times \vec A$ for some function $\vec A$. Plugging this into the third of Maxwell's equations gives
$$\nabla \times \vec E = -\frac{\partial}{\partial t}\left(\nabla \times \vec A\right)= -\nabla \times \left(\frac{\partial}{\partial t} \vec A\right)$$ so $$ \nabla \times \left(\vec E + \frac{\partial}{\partial t} \vec A\right) = 0$$ This implies that $$ \vec E + \frac{\partial}{\partial t}\vec A = -\nabla \phi$$ for some function $\phi$. From there, we find
$$\nabla^2 \phi + \frac{\partial}{\partial t}\left(\nabla \cdot \vec A\right) = -\rho/\epsilon_0$$ and $$\nabla \times \left(\nabla \times \vec A\right)=\nabla\left(\nabla \cdot \vec A\right) - \nabla^2 \vec A = \mu_0 \vec J + \epsilon_0 \mu_0 \left(-\frac{\partial}{\partial t}\left(\nabla \phi\right) - \frac{\partial^2}{\partial t^2} \vec A \right)$$
At this point, the answer is the same regardless of our choice of gauge, but the answer is most clear if we choose the Lorentz gauge condition:
$$\frac{\partial\phi}{\partial t} = -\frac{1}{\epsilon_0 \mu_0}\left(\nabla \cdot \vec A\right)$$
So the equations for the potentials become
$$\left(\frac{1}{c^2}\frac{\partial^2}{\partial t^2} - \nabla^2 \right)\phi = \rho/\epsilon_0$$
$$\left(\frac{1}{c^2}\frac{\partial^2}{\partial t^2} - \nabla^2 \right)\vec A = \mu_0 \vec J $$ where $c^2 = \frac{1}{\epsilon_0 \mu_0}$.
So we see that the potentials obey the wave equations whose sources are the charge density and current density, respectively.
You can make things even more concrete by transforming back to the fields explicitly. Because $$\vec E = -\nabla \phi - \frac{\partial }{\partial t} \vec A$$ it follows that the electric field obeys the wave equation:
$$\left(\frac{1}{c^2}\frac{\partial^2}{\partial t^2} -\nabla^2 \right)\vec E = -\left(\nabla \rho/\epsilon_0 + \mu_0\frac{\partial }{\partial t} \vec J\right)$$
This may seem complicated, but notice that the right hand side is just a function of space and time which is determined by the local presence and movement of electric charges.
From a mathematical standpoint, finding solutions to this equation is straightforward. If $$\left(\frac{1}{c^2}\frac{\partial^2}{\partial t^2} -\nabla^2 \right)\vec E = \vec f(\vec r, t)$$ for some function $\vec f$, and if we impose the boundary condition that $\vec E\rightarrow 0$ as $\vec r \rightarrow \infty$, then
$$\vec E(\vec r, t) = \int d^3r' \ \frac{\vec f\big(\vec r',t-\frac{|\vec r'-\vec r|}{c}\big)}{4\pi|\vec r' - \vec r|}$$
fits the bill quite nicely. Notice that if we want to find the field at position $\vec r$ and time $\vec t$, we look at the charges and currents at all other points $\vec r'$ as they were at the time $t_r = t-\frac{|\vec r'-\vec r|}{c}$. Therefore, there is a built-in time delay due to the fact that changes in the electric field propagate outward at speed $c$. That's why $t_r$ is called the retarded time, and this solution for $\vec E$ is called the retarded solution.
So, there you go. Isolated charges are not allowed to pop into existence, and changes in the electric field (and the magnetic field, by the same argument) propagate outward at $c$, so there is a time delay depending on where you are relative to the wiggling charges.
