  We proceed to obtain an exact analytical solution of the Brans-Dicke (BD)
equations for the spatially flat ($k=0$) Friedmann-Lamaitre-Robertson-Walker
(FLRW) cosmological model in both cases of the absence and presence of the
cosmological constant. The solution method that we use to solve the field
equations of the BD equations is called the "invariants of symmetry groups
method" (ISG-method). This method is based on the extended Prelle-Singer (PS)
method and it employs the Lie point symmetry, $\lambda$-symmetry, and Darboux
polynomials (DPs). Indeed, the ISG-method tries to provide two independent
first-order invariants associated to the one-parameter Lie groups of
transformations keeping ordinary differential equations (ODEs) invariant, as
solutions. It should be noted for integrable ODEs, the ISG-method guarantees
the extraction of these two invariants. In this work for the BD equations in
FLRW cosmological model, we find the Lie point symmetries, $\lambda$-symmetries
and DPs, and obtain the basic quantities of the extended PS method (which are
the null forms and the integrating factors). By making use of the extended PS
method we find two independent first-order invariants, in such a way
appropriate cosmological solutions from solving these invariants as a system of
algebraic equations are simultaneously obtained. These solutions are wealthy so
that they include many known special solutions, such as O'Hanlon-Tupper vacuum
solutions, Nariai's solutions, Brans-Dicke dust solutions, inflationary
solutions, and etc.
