Many-body eigenstates beyond the Gaussian approximation can be constructed in terms of local integrals of motion (IOMs), although their actual computation has been, until now, a daunting task. We present a practical computation of IOMs based on displacement transformations and apply it to the problem of many-body localization in disordered one-dimensional systems. It represents a general and systematic way to extend Hartree-Fock and configuration interaction theories to higher order. Our method combines minimization of energy and energy variance of a reference state with exact diagonalization. We show that our implementation for one-dimensional disordered fermions is able to perform ground-state calculations with high precision for relatively large systems. Since it keeps track of the IMOs forming a reference state, our method is particularly efficient dealing with excited states, both in accuracy and in the number of different states that can be constructed.