Recent advances in experimental physics opens up for improved time resolution on measurements of ionizing processes in atoms caused by an EM-field. There exist non-relativistic models that describe photo-ionization events involving two photons, for various atoms, where relativistic effects, e.g. spin-orbit interactions, are small. The ambition of this thesis project was to develop a computational software that includes relativistic effects by starting from the Dirac equation. Investigation of the photo-ionization events are done by calculating the two-photon matrix elements, taking into account essential many body effects using Random Phase Approximation with Exchange (RPAE) for the atomic response to the interaction with the first photon. During the development of the software it was discovered that the choice of dipole operator, length or velocity, had a large influence on the results when including the second photon, and it is concluded that the more straightforward calculations are done with the dipole operator in length gauge. The examined results demonstrate a program that works as desired; it agrees well with the non-relativistic calculations for the lighter elements and provides data for more detailed analysis of heavier elements.