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 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 for the first time 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.