  Recently considerable efforts have been devoted to computing cosmological
correlators and the corresponding wavefunction coefficients, as well as
understanding their analytical structures. In this note, we revisit the
computation of these ``cosmological amplitudes" associated with any tree or
loop graph for conformal scalars with time-dependent interactions in the
power-law FRW universe, directly in terms of iterated time integrals. We start
by decomposing any such cosmological amplitude (for loop graph, the
``integrand" prior to loop integrations) as a linear combination of {\it basic
time integrals}, one for each {\it directed graph}. We derive remarkably simple
first-order differential equations involving such time integrals with edges
``contracted" one at a time, which can be solved recursively and the solution
takes the form of Euler-Mellin integrals/generalized hypergeometric functions.
By combining such equations, we then derive a complete system of differential
equations for all time integrals needed for a given graph. Our method works for
any graph: for a tree graph with $n$ nodes, this system can be transformed into
the {\it canonical differential equations} of size $4^{n{-}1}$ quivalent to the
graphic rules derived recently%so-called ``kinematic flow", and we also derive
the system of differential equations for loop integrands {\it e.g.} of all-loop
two-site graphs and one-loop $n$-gon graphs. Finally, we show how the
differential equations truncate for the de Sitter (dS) case (in a way similar
to differential equations for Feynman integrals truncate for integer
dimensions), which immediately yields the complete symbol for the dS amplitude
with interesting structures {\it e.g.} for $n$-site chains and $n$-gon cases.
