The main problem in your proposed equation is that the electromagnetic equation with the D'Alambertian over the vector potential is not in Hamiltonian form, this means that the separation of solutions in Sturm–Liouville eigenstates of the energy operator is not manifest in the equation. Without that, you cannot find eigenstates of the coupled system.

You might find this dissertation interesting: On the canonical formulation of electrodynamics and wave dynamics

In there, the author analizes a Hamiltonian formulation for the electrodynamic field that is amenable to numerical solution coupled with the Schrödinger equations. Depending on what you actually want to find this should suffice (or not)

Regarding to your other question, i'm not confident giving an authoritative answer to that, so i'll let others jump on it

The main problem in your proposed equation is that the electromagnetic equation with the D'Alambertian over the vector potential is not in Hamiltonian form, this means that the separation of solutions in Sturm–Liouville eigenstates of the energy operator is not manifest in the equation. Without that, you cannot find eigenstates of the coupled system.

You might find this dissertation interesting: On the canonical formulation of electrodynamics and wave dynamics

In there, the author analizes a Hamiltonian formulation for the electrodynamic field that is amenable to numerical solution coupled with the Schrödinger equations. Depending on what you actually want to find this should suffice (or not)

Regarding to your other question, i'm not confident giving an authoritative answer to that, so i'll let others jump on it