  Viable models of modified gravity which satisfy both local as well as
cosmological tests are investigated. It is demonstrated that some versions of
such highly non-linear models exhibit multiply de Sitter universe solutions,
which often appear in pairs, being one of them stable and the other unstable.
It is explicitly shown that, for some values of the parameters, it is possible
to find several de Sitter spaces
  (as a rule, numerically); one of them may serve for the inflationary stage,
while the other can be used for the description of the dark energy epoch. The
numerical evolution of the effective equation of state parameter is also
presented, showing that these models can be considered as natural candidates
for the unification of early-time inflation with late-time acceleration through
dS critical points. Moreover, based on the de Sitter solutions, multiply SdS
universes are constructed which might also appear at the (pre-)inflationary
stage. Their thermodynamics are studied and free energies are compared.
