  A description of photoassociation by CW laser is formulated in the framework
of grid methods. The Hamiltonian describing one or several bound states coupled
to a multiple of continuum manifolds via a radiative field is written in the
energy representation and diagonalized. The generality of the treatment allows
to compute accurately and efficiently physical properties such as non-linear
high-intensity energy shifts, line shapes, and photoassociation rates both for
isolated and non-isolated resonances. Application is given to sodium
photoassociation in the experimental conditions of Mc Kenzie {\it et al} [Phys.
Rev. Lett. {\bf 88}, 090403 (2002)]. Inverted region for the dependency of the
rate vs. the intensity and non-symmetric lineshapes were predicted to occur
above the saturation limit. Comparison with the model of Bohn and Julienne
[Phys. Rev. A {\bf 60}, 414 (1999)] is discussed.
