  In high energy proton-nucleus collisions, the single- and double-inclusive
soft gluon productions at the leading order have been calculated and
phenomenologically studied in various approaches for many years. These studies
do not take into account the saturation and multiple rescatterings in the field
of the proton. The first saturation correction to these leading order results
(the terms that are enhanced by the combination $\alpha_s^2 \mu^2$, where
$\mu^2$ is the proton's color charge squared per unit transverse area) has not
been completely derived despite recent attempts using a diagrammatic approach.
This paper is the first in a series of papers towards analytically completing
the first saturation correction to physical observables in high energy
proton-nucleus collisions. Our approach is to analytically solve the classical
Yang-Mills equations in the dilute-dense regime using the Color Glass
Condensate effective theory and compute physical observables constructed from
classical gluon fields. In the current paper, the Yang-Mills equations are
solved perturbatively in the field of the dilute object (the proton).
Next-to-leading order and next-to-next-to-leading order analytic solutions are
explicitly constructed. A systematic way to obtain all higher order analytic
solutions is outlined
