The interplay between interactions and quenched disorder can result in rich dynamical quantum phenomena far from equilibrium, particularly when many-body localization prevents the system from full thermalization. With the aim of tackling this interesting regime, here we develop a semianalytical flow equation approach to study the time evolution of strongly disordered interacting quantum systems. We apply this technique to a prototype model of interacting spinless fermions in a random on-site potential in both one and two dimensions. Key results include (i) an explicit construction of the local integrals of motion that characterize the many-body localized phase in one dimension, ultimately connecting the microscopic model to phenomenological descriptions, (ii) calculation of these quantities in two dimensions, and (iii) an investigation of the real-time dynamics in the localized phase which reveals the crucial role of l-bit interactions for enhancing dephasing and relaxation.