  In this paper we develop further and refine the method of differential
equations for computing Feynman integrals. In particular, we show that an
additional iterative structure emerges for finite loop integrals. As a concrete
non-trivial example we study planar master integrals of light-by-light
scattering to three loops, and derive analytic results for all values of the
Mandelstam variables $s$ and $t$ and the mass $m$. We start with a recent
proposal for defining a basis of loop integrals having uniform transcendental
weight properties and use this approach to compute all planar two-loop master
integrals in dimensional regularization. We then show how this approach can be
further simplified when computing finite loop integrals. This allows us to
discuss precisely the subset of integrals that are relevant to the problem. We
find that this leads to a block triangular structure of the differential
equations, where the blocks correspond to integrals of different weight. We
explain how this block triangular form is found in an algorithmic way. Another
advantage of working in four dimensions is that integrals of different loop
orders are interconnected and can be seamlessly discussed within the same
formalism. We use this method to compute all finite master integrals needed up
to three loops. Finally, we remark that all integrals have simple Mandelstam
representations.
