At variance with fully inclusive quantities, which have been computed already at the two- or three-loop level, most exclusive observables are still known only at one-loop, as further progress was hampered so far by the greater computational problems encountered in the study of multi-leg amplitudes beyond one loop. We show in this paper how the use of tools already employed in inclusive calculations can be suitably extended to the computation of loop integrals appearing in the virtual corrections to exclusive observables, namely two-loop four-point functions with massless propagators and up to one off-shell leg. We find that multi-leg integrals, in addition to integration-by-parts identities, obey also identities resulting from Lorentz-invariance. The combined set of these identities can be used to reduce the large number of integrals appearing in an actual calculation to a small number of master integrals. We then write down explicitly the differential equations in the external invariants fulfilled by these master integrals, and point out that the equations can be used as an efficient method of evaluating the master integrals themselves. We outline strategies for the solution of the differential equations, and demonstrate the application of the method on several examples.