We outline a rigorous method which can be used to solve the many-body Schrödinger equation for a Coulomb interacting electronic system in an external classical magnetic field as well as a quantized electromagnetic field. Effects of the geometry of the electronic system as well as the polarization of the quantized electromagnetic field are explicitly taken into account. We accomplish this by performing repeated truncations of many-body spaces in order to keep the size of the many particle basis on a manageable level. The electron-electron and electron-photon interactions are treated in a nonperturbative manner using "exact numerical diagonalization." Our results demonstrate that including the diamagnetic term in the photon-electron interaction Hamiltonian drastically improves numerical convergence. Additionally, convergence with respect to the number of photon states in the joint photon-electron Fock space basis is fast. However, the convergence with respect to the number of electronic states is slow and is the main bottleneck in calculations.