Why does the Method of Successive Approximations for a Differential Equation work? Time dependent perturbation theory in quantum mechanics is often derived using the Method of Successive Approximations for a Differential Equation.  I have not seen an explanation or  a more rigorous outline of proof as to why this should work?
Normally we start with this equation:
$$
i\hbar \dot d_f=\sum_n\langle f^0|H^1(t)|n^0\rangle e^{i\omega_{fn} t}d_n(t)
$$
and then we say that to zeroth order the RHS is 0. and then we plug that into the first order and so on, as outlined in Shankar:

and as what $d_f(t)$ is. To zeroth order, we ignore the right-hand side of Eq. (18.2.5) completely, becasue of the explicit $H^1$, and get 
  $$\dot d_f=0 \tag{18.2.7}$$
  in accordance with our expectations. To first order, we use the zeroth-order $d_n$ in the right-hand side because $H_1$is itself of first order. This gives us the first-order equation
  $$
\dot d_f(t)=\frac{-i}\hbar\langle f^0|H^1(t)|i^0\rangle e^{i\omega_{fi} t} \tag{18.2.8}
$$
  the solution to which is, with the right initial conditions, is
  $$
d_f(t) = \delta_{fi} - \frac{i}\hbar\int_0^t \langle f^0|H^1(t')|i^0\rangle e^{i\omega_{fi} t'}\text dt' \tag{18.2.9}
$$

But how could we explain why this should work more rigorously?

For completeness I wanted to post an excerpt here from  Grynberg's Introduction to Quantum Optics which I though also answers the questions a bit more formally in terms of mathematics;

