Gfrktyl

Given a solution $A^{mu}(x)$ to Maxwell's equations
begin{equation}
Box A^{mu}(x)-partial^{mu}partial_{nu}A^{nu}=0tag{1}
end{equation}
which also satisfies some specified initial conditions at time $t_0$
begin{equation}
A^{mu}(vec{x},t_0)=f^{mu}(vec{x}),quad dot{A}^{mu}(vec{x},t_0)=g^{mu}(vec{x})tag{2}
end{equation}
we have that the function
begin{equation}
A^{'mu}(x)=A^{mu}(x)+partial^{mu}alpha(x)tag{3}
end{equation}
also satisfies the equations of motion, and if we arrange that the scalar function $alpha$ also satisfy that
begin{equation}
partial^{mu}alpha(vec{x},t_0)=0,quad partial^{mu}dot{alpha}(vec{x},t_0)=0 tag{4}
end{equation}
at the initial time $t_0$, then the new solution $A^{'mu}$ also satisfies the initial conditions. For example, the function
begin{equation}
alpha(vec{x},t)=(t-t_0)^5h(vec{x})e^{-(t-t_0)^2}
end{equation}
satisfies the conditions Eq.$(4)$ and also vanishes at $trightarrow pm infty$. Therefore, the solution to Eq.$(1)$ is not uniquely determined by the initial data Eq.$(2)$.

Question: If one simulates Eq.$(1)$ numerically on a computer, why is the field configuration at a later time not uniquely determined by the data in Eq.$(2)$?

Not all initial value problems have unique solution. Your example of $alpha$ function demonstrates that this initial value problem is of such kind.

In this case, the problem is in the system of partial differential equations

$$
partial_nupartial^nu A^{mu}-partial^{mu}partial_{nu}A^{nu}=0
$$
itself; it does not put enough constraint on the functions $varphi(mathbf x,t), mathbf A(mathbf x,t)$. It is somewhat similar to a situation in linear algebra that sometimes occurs where a system of $n$ linear equations for $n$ unknowns has infinity of solutions.

A slightly different way to see this: notice that nowhere in the above system of PDE can we find $partial_t^2 A^0$ or $partial_t A^0$ directly; only a spatial gradient of $partial_t A^0$ is present. The equations for $A^i$'s do not relate them directly with time derivatives of $varphi$.

This means that if we have a solution of the initial value problem $varphi(x,t),mathbf A(x,t)$ and replace the scalar potential by $varphi' = varphi(x,t)+ht^2$ at time $t = 0$ (where $h$ is a constant), the equations are still satisfied and at $t=0$, initial conditions are satisfied too. This would not be so obviously possible if the system contained directly time derivatives of $varphi$. Consider a slightly different system

$$
partial_nupartial^nu A^{mu}= 0,
$$
(which in EM theory can be derived as a result of the Lorenz gauge choice) - this does constrain $partial_t^2 varphi$, so the above argument fails. I think this system should have a unique solution, because it is very similar to a set of equations for independent harmonic oscillators. However, for proof better check with mathematicians.

Are you asking for the physical or mathematical explanation? Dan Yand's answer gives the physical explanation.

Regarding the mathematical question: On what basis would you expect the field configuration to be uniquely determined by its initial data? Unlike for (uncoupled) ODE's, there's no theorem to that effect for general linear homogeneous second-order PDEs.