This paper proposes a three-dimensional model for the evolution of micro and sub-micro scale flip chip solder joints. The coupled mechanisms of electromigration and the corresponding stress gradient are incorporated into a diffuse interface model. A semi-implicit Fourier spectral method and the preconditioned biconjugate-gradient method are proposed for the computation to achieve high efficiency and numerical stability. Our simulations demonstrate dynamic evolution of solder joints and breakages at the interface on the cathode side. It is also shown that the dynamically incorporated stress gradient considerably affects the evolution of solder joints counteracting the electromigration process.