numerical formulation of Dirac equation plus electromagnetic field I have the following equations describing the electron field in a (classic) electromagnetic field:
$$ c\left(\alpha _i\right.{\cdot (P - q(A + A_b) + \beta mc) \psi  = E \psi }  $$
where $A_b$ is the background field and $A$ is the one generated by the local Dirac field
I presume that the equation for the electromagnetic field $A$ generated by the electron would be:
$$\nabla_{\mu}\nabla^{\mu} A_{\nu} = \frac{\psi \gamma^{\nu} \psi}{\epsilon_0} $$
Question: Is there a way to numerically solve these systems of equations to find eigenstates of the system? 
Side Show Question: Are these eigenstates physically meaningful? do i still need to apply second quantization procedure in order to know which eigenstates are physically meaningful (i.e: stable) and which are not?
 A: 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
A: I guess these are tough questions, as the system is nonlinear. I can only give some references. In some of them, some numerical solutions of this system were found: Phys. Rev. A 60, 4291–4300 (1999) (also arXiv:physics/0001038 ), http://maths-old.anu.edu.au/research.publications/proceedings/039/CMAproc39-booth.pdf and references there. For what it's worth, in my work  arXiv:1111.4630 it is shown how to eliminate spinor field from the system (but complex electromagnetic potentials are introduced, which produce the same electromagnetic field, so their imaginary part is defined by one common function, so you have just 5 unknown real functions).
